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1.  SUMMARY 


A  mesoscale  numerical  modeling  system  is  adapted  for  use  into  the  lower 
stratosphere.  The  model  is  tested  and  an  automated  grid-nesting  scheme  is  developed  and 
tested  with  observed  initial  data.  The  purpose  of  the  model  and  the  new  nesting  scheme  is 
to  improve  the  prediction  of  regions  in  which  possible  severe  aviation  turbulence  occurs 
within  the  lower  stratosphere.  Simulation  results  show  a  possible  link  in  the  process  of 
stratospheric  turbulence  to  convectively  generated  large  amplitude  gravity  waves  that 
eventually  propagate  into  the  lower  stratosphere  and  break.  The  grid-nesting  scheme  is 
tested  to  determine  its  potential  utility. 

2.  INTRODUCTION 

The  operational  forecasting  of  tropospheric  turbulence  has  been  ongoing  for 
several  years.  Many  indices  are  generated  daily  from  operational  numerical  weather 
prediction  models.  The  National  Weather  Service  (NWS)  has  employed  the  Ellrod  Index 
(Ellrod  and  Knapp  1992),  the  National  Oceanic  and  Atmospheric  Administration 
Forecasting  Systems  Laboratory  (NOAA-FSL)  has  employed  indices  developed  by 
Marroquin  (1998)  based  on  turbulence  kinetic  energy  and  eddy  dissipation  rate,  and  the 
Research  Applications  Program  of  the  National  Center  for  Atmospheric  Research 
(NCAR-RAP)  employs  the  Graphical  Turbulence  Guidance  index  as  part  of  the  suite  of 
products  from  the  NWS  Rapid  Update  Cycle  (RUC)  II  model  (Sharman,  Wiener,  and 
Brown  2000).  These  systems  are  designed  for  use  in  the  troposphere  and  are  not  strongly 
coupled  to  the  model,  i.e.,  indices  are  simply  calculated  from  model-generated  dependent 
variables  and  any  grid  nesting  algorithms  are  independent  of  the  case  study  and  ongoing 
simulated  fields,  which  the  models  produce. 

The  processes  that  create  turbulence  in  the  stratosphere  are  likely  coupled  to  the 
troposphere  as  energy  regularly  is  exchanged  by  vertically  propagating  internal  gravity 
waves  between  the  two  atmospheric  regions.  This  vertical  coupling  is  greatly  facilitated 
during  moist  convection  or  the  development  of  large  amplitude  mountain  waves  as  these 
local  circulations  are  capable  of  modifying  their  immediate  environment  thus  resulting  in 
three-dimensional  fluxes  of  energy.  The  magnitude  of  the  vertical  flux  of  energy  is, 
therefore,  a  function  of  the  kinematic  and  thermodynamical  fields  that  evolve  in  response 
to  convective  motions,  i.e.,  the  perturbed  static  stability  and  vertical  wind  shear 
environments  surrounding  moist  or  dry  convection.  Simulating  stratospheric  turbulence  is 
an  extraordinarily  demanding  problem  because  of  this  vertical  coupling  as  well  as  the  fine 
scale  and  transient  nature  of  turbulence  itself.  Therefore,  to  develop  an  automated  system 
for  the  operational  prediction  of  lower  stratospheric  turbulence  potential  will  require  at 
least  three  fundamental  advances  in  mesoscale  numerical  weather  prediction:  1)  a  model 
with  significant  vertical  resolution  spanning  the  mid-upper  troposphere  and  lower 
stratosphere  in  order  to  resolve  the  vertical  coupling,  2)  a  model  with  a  “smart”  grid 
nesting  scheme  that  can  focus  on  the  regions  of  high  turbulence  potential  based  on  scale 
dependent  dynamical  indices  that  can  adapt  to  the  favorable  turbulence  forcing  processes 
that  vary  in  space  and  time  from  one  case  study  to  another,  and  3)  a  model  which  exploits 
as  much  of  the  observed  data  as  possible  in  its  initial  conditions  and  lateral  boundary 
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conditions  to  improve  the  accuracy  of  its  simulation.  We  are  modifying  a  modeling 
system  to  improve  the  automated  prediction  of  lower  stratospheric  regions  of  turbulence 
potential  in  which  we  incorporate  all  three  of  these  improvements.  During  year  one,  the 
focus  was  on  testing  the  quality  of  the  model  with  its  vertical  lid  extended  to  10  mb  and  its 
response  to  readily  available  conventional  gridded  data  sets  in  its  initial  state. 

Additionally,  work  commenced  on  formulating,  coding  and  testing  an  automated  grid 
nesting  scheme  or  autonest  based  on  an  original  algorithm. 

3.  METHODS,  ASSUMPTIONS  AND  PROCEDURES 

The  first  research  task  (Task  1)  involved  extending  the  model  lid  and  modifying 
the  vertical  structure  of  the  Non-Hydrostatic  Mesoscale  Atmospheric  Simulation  System 
(NHMASS)  model  version  6.3  (note  Table  1)  (e.g.,  Kaplan  et  al.  2000).  This  modeling 
system  was  turned  into  a  stratospheric  modeling  system,  i.e.,  STRATONHMASS  by 
extending  the  model  lid  to  10  mb,  increasing  the  number  of  vertical  levels  to  90  and 
modifying  its  preprocessor  to  ingest  observational  data  in  the  lower  stratosphere,  i.e.,  up  to 
~10  mb  or  ~30  km.  This  resulted  in  an  average  vertical  thickness  of  model  layers  of 
~300  m  with  a  nonuniform  vertical  focus  of  resolution  within  the  lower  stratosphere. 


2 


Table  1.  STRATONHMASS  (Version  6.3) 


iiiliiiiBiili 

ii 
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m 

i|g| 
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3-D  multivariate  01  procedure  used  to  blend  a  first  guess  fields  (e.g.,  previous 
MASS  simulation,  NCEP  model  output,  archived  Reanalysis  data)  with 
observations  from  a  variety  of  sensing  systems  (e.g.,  surface,  rawinsonde) 

Global  databases  of  terrain  height,  land  cover/land  use,  vegetation  index,  sea 
surface  temperature,  soil  texture,  snow  cover,  soil  moisture,  subsurface 
temperature,  and  sea  ice 


iSMlElll: 
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•  Option  of  hydrostatic  or  non-hydrostatic  primitive  equations  in  terrain-following  (op) 

vertical  coordinate  with  4* *-order  accurate  finite  differencing 

•  MPDATA  positive  definite  advection  scheme 

•  Option  of  one-way  or  two-way  interactive  nesting,  with  arbitrary  coarse/fine  grid 

spacing  ratio  and  unrestricted  number  of  nested  domains  for  one-way  nesting 


•  1.5  order  Turbulence  Kinetic  Energy  (TKE)  parameterization  with  surface  layer 

based  on  similarity  theory  formulation 

•  Surface  energy  budget  with  option  of  isothermal  or  non-isothermal  soil-vegetation 

canopy  formulations  and  heterogeneous  subgrid  scale  areas 

•  Surface  hydrology  includes  budget  equations  for  three  moisture  reservoirs  (cover 

layer,  shallow  and  deep  soil  layers)  and  snow  cover  which  incorporates  the 
effects  of  accumulation,  settling,  melting  and  sublimation 


•  Option  of  diagnostic  or  three  prognostic  schemes  with  varying  levels  of 
sophistication.  Prognostic  equations  for  cloud  water,  cloud  ice,  rain  water, 
snow  and  hail  using  a  bulk  microphysics  parameterization 


•  Option  of  Kuo-type  cumulus  parameterization  with  moist  downdraft  physics. 


Fritsch-Chappell  scheme,  Kain-Fritsch  scheme  or  Grell  scheme 


fKadiatibn 


si’ifiSiaslaigfSiii 


•  Longwave  and  shortwave  radiation  parameterized  in  surface  energy  budget  and 
in  the  free  atmosphere.  Interaction  with  atmospheric  water  vapor, 
liquid/frozen  water,  and  parameterized  sub-grid  clouds 
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The  second  research  task  (Task  2)  involved  testing  the  new  model  on  case  studies 
wherein  deep  moist  convection,  large  amplitude  gravity  waves  as  well  as  possible  terrain- 
induced  wave  phenomena  existed.  A  diverse  group  of  case  studies  was  simulated  in  an 
effort  to  test  the  reliability  of  the  model  with  its  vertically  extended  grid  as  well  as  the 
utility  of  conventional  gridded  and  rawinsonde  observations  in  the  lower  stratosphere. 

This  involved  producing  nested-grid  numerical  simulations  with  STRATONHMASS  for 
four  different  case  studies:  1)  a  convectively-generated  large  amplitude  gravity  wave 
propagating  on  a  well-developed  duct,  2)  a  winter  snowstorm  with  imbedded  gravity  and 
mountain  waves,  3)  a  tropical  system  with  possible  large  amplitude  gravity  waves  and 
deep  convection,  and  4)  a  winter  snowstorm  with  well-developed  maritime  convection 
and  possible  large  amplitude  gravity  waves.  These  case  studies  were  nested  employing 
the  newly  developed  autonesting  code  and  the  simulations  were  validated  against  all 
available  conventional  and  many  forms  of  asynoptic  data  including  wind  profilers, 

Doppler  radar  and  satellite  datasets  containing  Advanced  Microwave  Sounding  Unit 
(AMSU)  data.  Validation  indicated  reasonable  model  simulation  of  large-scale 
stratospheric  winds  and  temperatures  as  well  as  clouds  and  precipitation. 

Consistent  with  the  generation  of  the  four  sets  of  nested-grid  simulations,  the 
development  of  a  primitive  autonest  code  was  undertaken.  The  autonest  code 
development  was  accompanied  by  the  analyses  of  model-generated  dynamical  processes 
that  resulted  in  the  focusing  of  high  frequency  wave  energy  in  the  lower  stratosphere  both 
of  which  served  as  research  Tasks  3  and  4.  The  coarse  mesh  18-km  numerical  simulations 
were  nested  to  6-km,  2-km,  667-m,  222-m  and  71-m  horizontal  resolution.  The  first 
generation  of  the  autonest  code  employed  theNCSUl  index  (e.g.,  Kaplan  et  al.  2004). 

This  was  coded  to  nest  in  a  region  and  time  based  on  100-mb  NCSUl  threshold  values. 
The  second  generation  of  autonest  was  being  coded  to  nest  based  on  the  NCSUl  index  at 
all  grid  resolution  simulations  S2  km  and  then  to  employ  the  Scorer  Parameter  variation 
in  the  vertical  at  horizontal  resolutions  below  2  km.  The  third  generation  of  autonest  was 
being  coded  to  nest  based  on  threshold  values  of  eddy  dissipation  rate  (EDR)  from  the 
complete  turbulence  kinetic  energy  equation  at  all  grid  resolutions  below  222-m 
horizontal  resolution.  These  equations  are  listed  below. 


NCSUl  Index  =  [c/  •  V^l 

(1) 

U 

Scorer  Parameter  = — ^ - — 

(2) 

U 

EDR  = 


dx  dy  dz 


—r-,.du  d  .~i-s  1  d  .—i—f.  de 

-(uw  )-~  -^(w  e)-^~(wp  )-  — 
oz  oz  p  oz  ot 


(3) 
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Task  5  involved  initial  testing  of  code  that  will  be  used  to  improve  the  initial 
asynoptic  remotely-sensed  data  to  be  employed  in  the  model’s  initial  conditions.  This  will 
focus  on  relative  humidity  data  derived  from  infrared  cloud  top  temperatures. 

Task  6  involved  modifying  a  version  of  the  turbulence  kinetic  energy  equation  for 
use  in  the  autonesting  scheme  so  that  the  most  comprehensive  representation  of  the  eddy 
dissipation  rate  can  be  calculated.  This  involved  adding  the  advective  terms  to  make  the 
equation  complete,  as  is  described  in  Appendix  B. 

Finally,  Task  7  involved  several  required  Air  Force  Research  Laboratory  reports, 
two  conference  preprints  and  two  journal  articles,  which  were  prepared  or  were  in  the 
process  of  being  prepared  at  the  end  of  year  one  in  an  effort  to  document  the  research. 

4.  RESULTS  AND  DISCUSSION 

As  noted  in  Appendices  A  and  B,  the  paradigm  developed  to  modify  and 
implement  successive  generations  of  the  autonest  code  is  based  on  simulated  dynamical 
adjustment  processes  in  the  lower  stratosphere.  Two  case  studies,  in  particular,  have  been 
used  to  test  the  paradigm,  i.e.,  the  12-13  December  2002  convective  gravity  wave  event 
and  the  16-17  February  2003  Mid- Atlantic  extreme  snowfall  event.  Additionally,  the 
landfall  of  tropical  storm  Gaston  during  30-3 1  August  2004  and  the  24-25  January  2000 
Mid-Atlantic  coastal  storm  were  also  being  simulated.  In  both  extensively  simulated  case 
studies,  i.e.,  12-13  December  2002  and  16-17  February  2003,  large  amplitude  internal 
gravity  waves  occur,  due  in  part  to  the  development  of  an  unbalanced  jet  streak  exit  region 
(e.g.,  Uccellini  and  Koch  1987;  Kaplan  et  al.  1997). 

The  16-17  February  2003  snowstorm,  and  its  potential  implications  for 
stratospheric  turbulence,  is  discussed  in  detail  in  Appendix  A.  In  this  section  of  the  report 
we  will  describe  the  large  amplitude  gravity  wave  case  study  of  12-13  December  2002 
wherein  mesoscale  convective  systems  developed  above  a  stable  layer.  The  turbulence 
kinetic  energy  diagnostics  calculated  for  this  case  study  are  depicted  in  Appendix  B.  In 
this  section  we  will  focus  on  the  scale  contraction  process  that  results  in  stratospheric 
wave  amplification  and  breaking,  necessitating  the  turbulence  kinetic  energy  analysis  in 
Appendix  B. 

Figure  1  depicts  the  radar  imagery  valid  at  2315  UTC,  12  December  2002.  The 
large  convective  complex  system  located  across  eastern  Texas  and  Louisiana  originated 
over  the  southern  Texas  Gulf  coast  at  about  1630  UTC,  12  December.  Collocated  with 
this  convective  complex  and  its  propagation  across  eastern  Texas  north  of  a  stable  layer,  is 
the  large  amplitude  gravity  wave  diagnosed  in  the  NOAA  Global  Positioning  System 
(GPS)  microbarogram  in  Figure  2  at  Palestine  (PAL),  in  southeastern  Texas.  The 
microbarogram  indicates  a  wave  of  depression  associated  with  a  nearly  6-mb  pressure  fall 
in  about  two  hours  just  before  2230  UTC.  This  fall  is  a  likely  signal  of  a  large  amplitude 
gravity  wave  based  on  this  single  station  trace  as  well  as  other  GPS  traces  in  eastern  Texas 
and  Louisiana,  which  show  similar  temporally  coherent  signals.  Thus,  even  in  the 
absence  of  a  detailed  surface  covariance  analysis,  inferential  observational  proof  of  the 
wave  exists  which  will  be  bolstered  by  the  mesoscale  simulations. 
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Figure  1.  NWS  radar  summary  valid  at  2315  UTC,  12  December  2002. 


Figure  2.  GPS  mean  sea  level  pressure  (mb)  at  Palestine,  Texas,  valid  from  12/0100  UTC 
-  15/0000  UTC,  December  2002. 


Winds  from  profiling  radars  at  Palestine  are  depicted  in  Figures  3  (a)  and  (b)  at 
hourly  and  6-minute  time  intervals.  The  time  of  wave  passage  in  the  pressure  data  was 
shortly  before  2230  UTC.  It  is  during  this  period  that  the  6-  to  1 0-km  winds  in  the  hourly 
data  show  acceleration  and  turning  from  southwest  to  south.  Above  this  level,  there  is  a 
brief  acceleration  to  the  west,  particularly  at  about  150  mb  or  ~1 4  km,  which  is  in  the 
lower  stratosphere  immediately  above  the  tropopause.  The  6-minute  data  confirms  an 
increase  from  the  west  and  then  decrease  of  westerly  fiow  at  that  level  after  about  2330 
UTC  when  the  wave  of  depression  was  located  well  east  of  Palestine  in  microbarograms 
in  southwestern  Louisiana. 
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Figure  3(a).  Winds  (m/s)  from  vertically  profiling  radars  at  PalestinCj  Texas,  every  hour 
1600  UTC  -  0300  UTC,  12-13  December  2002,  from  right  to  left. 


Figure  3(b).  Winds  (m/s)  from  vertically  profiling  radars  at  Palestine,  Texas,  every  6 
minutes  from  2254  UTC  -  0000  UTC,  12-13  December  2002,  from  right  to  left. 
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The  inferences  to  be  drawn  from  all  these  data  are:  1)  that  convection  was  fairly 
closely  related  to  the  wave  of  depression  and  2)  that  the  wave/convective  system  likely 
played  some  role  in  modifying  the  winds  aloft  above  and  below  the  tropopause.  This  can 
be  seen  in  the  accelerating  of  the  flow  below  the  tropopause,  followed  by  a  deceleration  in 
velocities  below  the  tropopause,  as  well  as  accelerations  above  the  tropopause  nearly 
coincident  with  the  lower  layer  deceleration.  This  is  analogous  to  the  wave  and  its 
attendant  convective  heating  acting  to  vertically  expand  the  height  of  the  previous 
background  jet  stream  flow  or  to  force  a  vertical  mass  convergence  between  6  and  14  km. 
This  process,  therefore,  should  likely  result  in  the  evacuation/decrease  of  horizontal 
momentum  at  300  mb  and  increase  of  momentum  at  150  mb  on  or  about  0000  UTC,  13 
December  2002.  The  consequence  is  that  southerly  momentum  is  being  forced  upwards 
from  the  lower  troposphere  and  westerly  momentum  is  being  forced  upwards  from  the 
upper  troposphere  being  forced  into  the  lower  stratosphere  while  the  down  and  upstream 
flows  are  being  modified  as  well.  This  restructuring  of  the  flow  can  be  further  confirmed 
by  viewing  rawinsonde-based  gridded  wind  analyses  depicted  in  Figures  5  and  6.  Thus, 
the  wave  and  its  attendant/nearly  collocated  diabatic  forcing  are  acting  to  vertically 
expand  the  level  of  momentum  shears  both  upward  and  outward  in  3-dimensional  space 
from  their  pre-convective  levels  in  the  vertical  and,  therefore,  acting  to  modify  the  wind 
shears  in  the  lower  stratosphere.  This  is  all  accompanying  the  convectively  induced 
energy  fluxes  that  were  hypothesized  earlier  to  be  critical  to  generating  stratospheric  wave 
breaking. 

In  an  effort  to  test  the  new  stratospheric  model  and  the  first  generation  autonest 
code,  a  sequence  of  numerical  simulations  were  performed  with  the  model  employing  a 
three-dimensional  grid  matrix  of  162x162x90  points.  Figure  4  depicts  the  6-km,  2-km, 
667-m,  222-m,  and  71-m  regions  of  model  integration  for  this  12-13  December  2002  case 
study,  all  of  which  were  embedded  within  the  18-km  grid,  which  covered  much  of  the 
continental  U.S.  The  model  coarse  mesh  grid  was  initialized  at  0000  UTC,  12  December 
2002,  from  NCEP  Reanalysis  gridded  data  and  reanalyzed  rawinsonde  and  surface 
observations.  Each  finer  mesh  simulation  was  integrated  for  a  progressively  shorter 
period  of  time  centered  on  and  consistent  with  wave  passage  near  the  Texas/Louisiana 
border  region.  The  initial  conditions  and  lateral  boundary  conditions  were  derived  from 
each  successively  coarser  mesh  model  simulation  for  a  one-way  nested-grid  configuration. 
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Figure  4.  Grid  mesh  locations  for  the  6-km,  2-km,  667-m,  222-m  and  71 -m 
STRATONHMASS  model  runs. 


Figures  5  and  6  depict  the  18-km  wind  validation  against  the  rawinsonde 
observations  at  300  and  150  mb.  Most  notable  in  the  simulations  is  the  previously- 
referenced  decelerating  flow  at  300  mb  and  accelerating  flow  at  150  mb  near  the 
Texas/Louisiana  border  at  0000  UTC,  13  December.  This  is  deduced  from  the  apparent 
300-mb  wind  minimum  near  the  convection  in  east  Texas  and  maximum  in  the  jet 
entrance  region  aloft  at  1 50  mb  in  the  same  region.  The  timing  is  consistent  with  gravity 
wave  and  mesoscale  convective  system  passage  through  this  region  Just  southeast  of 
Palestine,  Texas.  The  model  replicates  the  general  pattern  of  stronger  winds  at  150  mb 
relative  to  that  at  300  mb  consistent  with  observed  data  in  this  region.  As  such,  the 
vertical  wind  shear  above  the  tropo pause,  i.e.,  within  the  25-  to  50-mb  deep  layer  centered 
on  150  mb,  is  likely  very  large  due  to  the  expansion/vertical  and  horizontal  convergence 
of  the  momentum  field  above  the  convective  heating.  This  increased  wind  shear  above 
the  wave  is  consistent  with  the  6-km  simulated  300-mb  winds  in  Figure  7  wherein  the 
momentum  has  been  depleted  at  300  mb  and  forced  upward  and  outward  into  a  new 
highly  ageostrophic  jet  entrance  region  to  the  north  and  west  of,  as  well  as  above,  the 
simulated  and  observed  convection.  Hence,  the  diabatic  heating  is  forcing  kinetic  energy 
fluxes  both  into  the  stratosphere  and  upstream  to  the  north  and  west  of  the  surface  wave 
signal  in  the  simulated  fields.  This  can  be  seen  by  comparing  Figures  7  and  8  with  the  2- 
km  nested-grid  mean  sea  level  pressure  trough  and  ridge  accompanying  the  wave  near  the 
Texas/Louisiana  border  region  and  its  highly  ageostrophic  surface  flow,  which  lies  just 
southeast  of  the  depleted  flow  at  300  mb.  Hence,  tilting  backward  from  the  simulated 
surface  wave  is  a  kinetic  energy  minimum  at  300  mb  and  maximum  above  the  tropopause 
reflecting  the  phase  shift  of  kinetic  energy  adjustments  to  the  convectively-induced 
heating. 
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Figure  5.  Comparison  of  1 8-km  STRATONHMASS  simulated  300-mb  isoheights  (m) 
and  isotachs  (knots)  [top]  with  analyzed  observed  isotachs  (knots)  [bottom]  valid  at  0000 
UTC,  1 3  December  2002. 
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Figure  6.  Comparison  of  18-km  STRATONHMASS  simulated  1 50-mb  isoheights  (m) 
and  isotachs  (knots)  [top]  with  analyzed  observed  isotachs  (knots)  [bottom]  valid  at  0000 
UTC,  13  December  2002. 
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Figure  7.  6-km  STRATONHMASS  300-mb  isoheights  (m)  and  isotachs  (knots)  valid  at 
2300  UTC,  12  December  2002. 
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Figure  8.  2-km  STRATONHMASS  simulated  mean  sea  level  pressure  (mb)  and  surface 
winds  (knots)  valid  at  2300  UTC,  12  December  2002. 


Further  simulated  proof  of  this  sloping  shear  zone  and  its  implication  for 
stratospheric  turbulence  can  be  gleaned  from  Figures  9  and  10,  which  depict  very  high 
resolution,  i.e.,  222  m,  simulated  precipitation  and  surface  winds  as  well  as  the  100  mb 
NCSUl  index.  Note  how  the  turbulence  index  at  100  mb  is  displaced  well 
upstream/northwest  from  the  precipitation.  This  signal  of  large  index  values  is  consistent 
with  the  discontinuous  advection  of  stream  wise  momentum,  vorticity  gradients,  vertical 
wind  shear  and  decreased  static  stability  maximizing  many  kilometers  to  the  northwest  of 
the  border  and  to  the  northwest  of  the  ongoing  simulated  precipitation.  Thus,  the  region 
of  potential  horizontal  flow  bifurcation  and  vertically  propagating  wave  structures  may 
maximize  well  displaced  in  space  and  time  from  the  actual  convective  heating.  The 
horizontal  flow  bifurcation  and  subsequent  vertically  propagating  gravity  wave  energy 
resulting  from  the  shear  and  instability  are  likely  caused  by  the  temporally  integrated 
effects  of  convection  on  the  environmental  dependent  variables  rather  than  at  the  location 
of  the  convection  itself. 


Figure  9.  222-m  STRATONHMASS  simulated  1 -minute  precipitation  (mm)  and  surface 
winds  (knots)  valid  at  2230  UTC,  12  December  2002. 


Figure  10.  222-m  STRATONHMASS  simulated  100-mb  NCSUl  index  (s'^  x  10'®)  valid 
at  2230  UTC,  12  December  2002. 
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To  visualize  the  consequence  of  these  processes  better  at  the  finest  scales  of 
motions  a  vertical  cross  section  corresponding  to  cross  section  M  depicted  in  Appendix  B, 
is  oriented  west  to  east  as  can  be  seen  with  simulated  isentropes  and  isotachs  depicted  in 
Figures  11-13  for  three  time  periods  separated  by  three  minutes  each  from  0015  -0021 
UTC.  The  location  of  the  cross  section  bisects  the  high  values  ofNCSUl  depicted  in 
Figure  10  along  a  west  -  east  axis.  The  isentropic  perturbation  above  the  wind  shear  zone 
near  90  mb  is  consistent  with  a  vertically  varying  value  of  the  Brunt-Vaisalla  frequency 
and  the  second  derivative  of  wind  with  respect  to  height  just  above  the  tropopause.  These 
shear  and  stability  terms  control  vertical  fluctuations  in  the  stability  and  wind  shear 
consistent  with  a  transition  region  between  wave  trapping  and  critical  layer  formation,  i.e., 
where  the  value  of  the  Scorer  parameter  changes  signs  from  positive  to  negative  in  the 
vertical.  Critical  layer  theory  dictates  that  wave  breaking  would  be  ongoing  where  the 
Scorer  parameter  decreases  to  less  than  zero  in  such  a  region.  The  wave  amplification  and 
overturning  seen  in  Figures  11-13  are  consistent  with  shearing  instability  occurring  where 
a  convective  outflow  jet  has  modified  the  lower  stratosphere.  The  model  is  predicting 
wave  breaking  here,  hence  it  is  a  region  where  the  autonesting  code  should  sense  the 
increase  in  Scorer  parameter  variation  and  its  possible  organization  of  local  shear  and 
buoyancy  maxima  ideal  for  the  generation  of  turbulence  kinetic  energy.  The 
correspondence  between  this  vertically  varying  zone  of  Scorer  parameter  and  increasing 
eddy  dissipation  rate  can  be  seen  depicted  at  the  grid  point  in  Figure  14.  At  approximately 
90  mb,  i.e.,  the  level  of  wave  breaking,  the  local  increase  in  turbulence  kinetic  energy  wilt 
force  the  model  autonest  code  to  flag  this  location  as  ideal  at  this  time  for  nesting  to  71  m. 


Figure  11.  222-m  STRATONHMASS  vertical  cross  section  from  30.42  N,  -93.92  W 
and  30.42  N,  -93.67  W.  Isentropes  (K)  and  isotachs  (knots)  valid  at  2224  UTC, 

12  December  2002. 
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Figure  12.  Same  as  Figure  1 1,  but  at  2227  UTC,  12  December  2002. 
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Figure  14.  222-m  STIIATONHMASS  Scorer  parameter  (solid;  m'^)  and  eddy  dissipation 
rate  (dashed;  m^s'^)  valid  at  the  middle  point  of  cross  section  in  Figure  13  at  2233  UTC, 
12  December  2002. 


The  nested  grid  71-m  cross  section  (M2)  of  simulated  isentropes  and  isotachs 
oriented  nearly  orthogonal  to  cross  section  M  in  Figures  1 1-13  is  depicted  in  Figure  15  at 
0024  and  0025  UTC,  12  December  2002.  The  maximum  vertical  variation  of  Scorer 
parameter  in  the  stratosphere  near  90  mb  can  be  seen  in  Figure  16  with  the  superimposed 
vertical  variation  of  eddy  dissipation  rate  at  the  same  time  and  level  from  the  71-m 
simulation.  This  orthogonal  cross  section  depicts  shallower  cross-wave  instabilities  along 
the  periphery  of  the  feature  in  Figures  11-13.  Such  possible  regions  of  wave  breaking  and 
turbulence  may  be  validated  against  AMSU  data  in  the  near  future  if  such  data  is  available 
from  NOAA. 


17 


Figure  15.  Same  as  Figure  1 1 ,  but  at  71 -m  resolution  orthogonal  to  cross  section  in 
Figure  13,  at  2238  UTC  [top]  and  2239  UTC  [bottom],  12  December  2002. 
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Figure  16.  Same  as  Figure  14,  but  from  the  71-m  simulation  valid  at  2239  UTC, 
12  December  2002. 


5.  CONCLUSIONS 

The  first  year  of  research  sho\ved  major  progress  towards  the  development  of  a 
stratospheric  real-time  turbulence  modeling  system.  During  this  period,  the  NHMASS 
model  was  tested  with  the  lid  raised  well  into  the  stratosphere  on  four  research  case 
studies.  Results  of  the  simulated  fields  were  satisfactory  when  validated  against 
conventional  and  asynoptic  observations.  An  autonest  algorithm  was  formulated  and 
tested  on  these  four  case  studies.  The  abili  ty  of  the  model  to  simulate  fine  scale 
convective  features  and  wave  phenomena,  as  well  as  the  autonest  to  direct  the  model  grid 
nests  to  finer  scales  in  regions  of  convectively  modified  circulations,  indicates  the 
potential  for  use  in  turbulence  generating  regions.  However,  validation  of  the  rhodeling 
system  lags  its  development  and  testing  due  to  the  lack  of  stratospheric  turbulence  ' 
observations. 

6,  RECOMMENDATIONS 

There  is  a  need  for  information  concerning  observed  turbulence  reports  in  the 
lower  stratosphere  in  the  existing  simulated  and  future  simulated  case  studies  to  validate 
the  modified  and  evolving  modeling  system. 
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Appendix  A 

Inertia-Gravity  Wave  Breaking  and  Lower-Stratospheric  Turbulence  During  the 

2003  Presidents’  Day  Event 


1.  INTRODUCTION 

Atmospheric  turbulence  is  both  a  significant  threat  to  aviation  throughout  the 
United  States  and  an  extremely  difficult  phenomenon  to  observe  and  forecast,  making 
turbulence  one  of  the  more  critical  research  topics  in  atmospheric  science.  Between  1986 
and  1993,  close  to  1000  people  were  injured  as  a  result  of  turbulence  impacting  aircraft 
(Cowen  1998),  and  for  commercial  airliners  alone,  tens  of  passengers  are  injured  each 
year  (Lane  et  al.  2003).  The  small  spatial  and  temporal  scale  of  turbulence  limits  both  the 
observation  and  accurate  forecasting  of  atmospheric  turbulence.  Cetola  (2003)  points  out 
that  the  identification  of  large-scale  precursors  of  turbulence  is  a  necessary  step  in  order  to 
predict  turbulence  using  the  relatively  course  observational  data  commonly  available. 
Turbulence  has  been  shown  to  have  a  wide  variety  of  sources,  including  mountain  wave 
breaking,  wind  shear-induced  Kelvin-Helmholtz  instabilities,  overshooting  convective 
updrafts,  and  inertia-gravity  wave  breaking.  The  purpose  of  this  study  is  to  relate  the 
development  of  inertia-gravity  waves  to  the  simulated  development  of  stratospheric 
turbulence  over  portions  of  West  Virginia  and  western  Maryland  during  the  2003 
Presidents’  Day  winter  storm  of  16-17  February  2003.  Since  limited  observational  data  is 
available  to  validate  the  existence  of  apparent  inertia-gravity  wave  activity  over  the  study 
region,  the  majority  of  this  study  makes  use  of  a  mesoscale  numerical  model,  the  Non- 
Hydrostatic  Mesoscale  Atmospheric  Simulation  System  (NHMASS)  version  6.3  (Kaplan 
et  al.  2000)  to  evaluate  the  complex  non-linear  dynamics  resulting  in  the  transition  of 
mesoscale  wave  activity  to  turbulence. 

Inertia-gravity  wave  sources  have  been  shown  to  include  shearing  instability, 
convection,  and  geostrophic  adjustment  (Uccellini  and  Koch  1987).  The  latter  mechanism 
will  be  evaluated  in  this  paper  due  to  the  elimination  of  the  first  two  mechanisms  in 
analyses  of  NHMASS  fields  (not  shown),  and  the  presence  of  a  strongly  unbalanced 
subtropical  jet  (STJ)  exit  region  immediately  upstream  of  the  region  of  propagating  wave 
activity.  Methods  of  distinguishing  inertia-gravity  wave  activity  from  other  mesoscale 
banded  phenomena,  including  quadrature  between  mid-tropospheric  potential  temperature 
(0)  and  vertical  velocity  (w)  fields,  and  evaluating  the  in-phase  relationship  between  the 
surface  fields  of  perturbation  pressure  (p’)  and  wave-normal  horizontal  wind  speed  (u*) 
(Koch  and  Golus  1988)  will  be  evaluated  to  partially  validate  the  classification  of  the 
phenomena  observed  as  inertia-gravity  wave  activity. 

High-resolution  simulations  will  be  used  to  analyze  the  transition  to  turbulence  in 
the  lower  stratosphere,  evaluating  the  process  of  wave  overturning,  generation  of  fine- 
scale  instabilities,  and  the  rapid  development  of  turbulence.  Lastly,  a  number  of 
turbulence  measurements,  including  turbulence  kinetic  energy  (TKE)  dissipation  and  the 
NCSUl  index,  an  index  sensitive  to  horizontal  gradients  of  momentum  and  low 
Richardson  number  values,  will  be  utilized  in  assessing  turbulence  in  the  simulations. 
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2.  MODEL  DESCRIPTION 


Due  to  the  lack  of  observational  data  available  at  a  scale  necessary  to  resolve  the 
mechanisms  discussed  in  this  paper,  a  mesoscale  model  is  employed  for  dynamical 
analysis.  The  Non-Hydrostatic  Mesoscale  Atmospheric  Simulation  System  (NHMASS) 
model  version  6.3  (Kaplan  et  al.  2000),  adapted  for  modeling  of  stratospheric  dynamics,  is 
used  for  this  study.  TKE  PBL  physics  (Therry  and  Lacarrerre  1983),  mixed-phase 
moisture  physics  (Lin  et  al.  1983;  Rutledge  and  Hobbs  1983),  and  the  Kain-Fritsch 
cumulus  parameterization  scheme  (Kain  and  Fritsch  1993)  were  used  in  all  simulations 
performed.  A  nested-grid  approach  was  used,  with  the  initial  grid  featuring  I8-km  grid 
spacing  and  each  subsequent  nest  finer  in  scale  by  a  factor  of  three.  The  fmest-scale 
horizontal  grid  spacing  used  was  222  m,  which  will  be  utilized  in  evaluating  stratospheric 
wave  breaking  and  generation  of  turbulence.  A  summary  of  simulations  performed  is 
found  in  Table  Al .  See  Figure  Al  for  a  schematic  of  grid  meshes  discussed  in  this  paper. 

3.  ANALYSIS  OF  FLOW  IMBALANCE  AND  IGW  GENERATION 

To  understand  the  unbalanced  flow  inertia-gravity  wave  energy  source,  one  needs 
to  first  consider  the  flow  characteristics  of  a  balanced  jet  streak,  depicted  in  a  schematic 
given  in  Figure  2.  One  notes  an  ageostrophic  component  directed  toward  lower 
streamfunction  (akin  to  geopotential  height  on  an  isobaric  surface)  in  the  entrance  region 
of  the  jet  streak,  and  a  component  directed  toward  higher  streamfunction  in  the  exit 
region.  The  direction  of  these  ageostrophic  components  is  consistent  with  acceleration 
(deceleration)  of  the  total  wind  in  the  entrance  (exit)  region  of  the  geostrophic  wind 
maximum,  and  represents  an  attempt  by  the  atmosphere  to  maintain  mass-momentum  (i.e., 
thermal  wind)  balance.  Any  deviation  from  this  simple  conceptual  model,  which  could 
disrupt  thermal  wind  balance,  will  be  considered  unbalanced flow,  a  classification  of  great 
importance  since  such  flow  is  known  to  generate  unbalanced  mesoscale  circulations,  such 
as  inertia-gravity  waves. 
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Table  Al.  Summary  of  MASS  Simulations  Performed 


Initialized 

(UTC) 

Resolution 

(km) 

^^99 

Modifications 

Hydrostatic 
/  Non- 
Hydrostatic 

Simulation 

Name 

02/14/03 

1200 

36 

72 

None 

H 

FULL36 

02/14/03 

1200 

36 

95,114, 90 

72 

Dry* 

H 

02/15/03 

1200 

18 

162,162,90 

60 

None 

NH 

02/15/03 

1200 

18 

162,162,90 

60 

Smooth  Terrain 

NH 

SMOOTH  18 

02/16/03 

1800 

6 

18 

None 

NH 

02/16/03 

1800 

6 

153,100,90 

18 

Smooth  Terrain 

NH 

SMOOTH6 

02/l'6/03 

1800 

6 

153,100,90 

18 

Dry 

NH 

02/16/03 

1830 

2 

200,125,90 

9 

None 

NH 

02/16/03 

1830 

2 

200,125,90 

9 

Smooth  Terrain 

NH 

SMOOTH2 

02/17/03 

0100 

0.670 

127,127,90 

2 

None 

NH 

FULL670 

02/17/03 

0100 

0.670 

127,127,90 

2 

Smooth  Terrain 

NH 

SMOOTH670 

02/17/03 

0100 

0.222 

162,162,90 

0.75 

None 

NH 

FULL222 

*  DRY  runs  differ  from  FULL  simulations  in  that  diabatic  heat  release  due  to  evaporation 
and  condensation  has  been  omitted.  Sensible  heating  remains  intact  in  all  simulations. 


b 


Figure  A1 .  (a)  36-km  Domain  (light  shading),  18-km  and  6-lcm  grid  meshes  (medium  and 
dark  shading,  respectively),  and  (b)  6-km  Domain  (darkest  shading),  2-km,  670-m,  and 
222-m  grid  meshes  (lightest  shading). 
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Figure  A2.  Schematic  representation  of  the  ageostrophic  motions  (heavy  arrows)  and 
associated  patterns  of  convergence  (CON)  and  divergence  (DIV)  in  the  vicinity  of  a 
straight  jet  streak  in  the  absence  of  along-contour  thermal  advection.  Assumed  to  be  near 
level  of  maximum  wind,  where  the  horizontal  wind  distribution  is  most  distinct  and  the 
flow  is  approximately  horizontal.  Solid  lines  indicate  geopotential  height  of  a  constant 
pressure  surface;  dashed  lines  are  isotachs  with  maximum  wind  speed  shaded.  [From 
Keyser  and  Shapiro  (1981)  and  Keyser  and  Shapiro  (1986)]. 

Figure  A3  is  a  300-hPa  analysis  of  U.S.  rawinsonde  data  valid  1200  UTC,  16  Feb 
2003,  and  indicates  a  component  of  the  total  wind  directed  toward  lower  streamfunction  in 
the  exit  region  of  the  (geostrophic)  STJ  streak  (indicated  by  the  oval  in  Figure  A3).  This 
is  clearly  inconsistent  with  the  balanced  jet  schematic  in  Figure  A1 .  Evidence  of 
unbalanced  flow  within  the  NHMASS  simulations  performed  is  given  in  Figure  A4,  an 
analysis  of  the  Lagrangian  Rossby  number  (Rol)  within  the  36-km  full  physics  NHMASS 
simulation  (FULL36)  at  (a)  16/1200,  (b)  16/1800,  and  (c)  17/0000,  where 

Rol  =  |DV/dt|  /  |fV|.  (Al) 

Caution,  though,  must  be  exercised  in  the  use  of  Rol  in  assessing  imbalance  in  highly 
curved  flow  such  as  that  seen  in  Figure  A3,  as  Rol  includes  the  effects  of  curved  flow, 
which  itself  can  be  described  through  gradient  wind  balance.  For  this  reason,  precise 
magnitudes  of  Rol  should  not  be  gleaned  from  Figure  A4;  only  an  evaluation  of  regions 
of  likely  flow  imbalance  should  be  made.  Large  values  of  Rol(>  0.5)  have  been  shown  to 
exist  for  flows  that  cannot  be  described  by  either  geostrophic  or  gradient  wind  balances, 
i.e.,  unbalanced  flow  (Uccellini  and  Koch  1987;  Koch  and  Dorian  1988;  Rozumalski  1998 
among  others),  and  will  be  utilized  here  to  confirm  the  unbalanced  flow  signal  in  the 
simulation.  One  notes  that  the  region  of  high  Ro  (~0.7-0.8)  at  16/1200  [Figure  A4(a)] 
over  much  of  Kentucky  and  southern  Indiana  is  part  of  the  same  region  inferred  as 
unbalanced  in  the  300-hPa  upper-air  analysis  (Figure  A3).  By  17/0000  [Figure  A4(c)],  the 
region  of  large  Ro  has  expanded  and  moved  to  the  northeast,  with  the  largest  Ro  (>0.8) 
over  northern  West  Virginia.  The  movement  of  this  region  of  apparently  unbalanced  flow 
is  consistent  with  the  motion  of  the  STJ  exit  region  in  both  the  upper-air  analysis  and 
simulated  300-hPa  fields  (not  shown). 
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Figure  A3.  300-hPa  NMC  analysis  valid  0000  UTC,  2/17/03.  Oval  denotes  region 
described  in  text. 


26 


Figure  A4.  FULL36  332  K  Lagrangian  Rossby  valid  (a)  16/1200,  (b)  16/1800,  and 
(c)  17/  0000.  The  star  in  (c)  marks  location  of  western  Maryland. 


With  the  presence  of  unbalanced  flow  in  the  STJ  exit  region  now  established  in 
both  the  observational  and  NHMASS  datasets,  an  assessment  of  IGW  activity  will  be 
made.  In  Figure  A5,  which  depicts  an  NOAA  Global  Positioning  System  (GPS)  surface 
pressure  trace  from  the  Louisville,  Kentucky  (LOU)  GPS  site,  one  notices  a  roughly  one- 
hPa  wave  of  elevation  near  16/2100  followed  one  hour  later  by  a  wave  of  depression  of 
opposite  amplitude.  This  pair  of  high-frequency  waves  is  embedded  within  a  much  longer 
period  (~48hr)  oscillation  associated  with  the  passage  of  surface  low  pressure. 

Reasonably  consistent  with  this  feature  in  the  GPS  pressure  trace  is  a  wave  of  depression 
noted  within  the  FULL  18  surface  pressure  tendency  fields  at  16/2000  from  central 
Kentucky  to  upstate  South  Carolina  [Figure  A6(a)].  It  should  be  noted  that  estimations  of 
horizontal  wavelength  and  vertical  wavelengths  of  94  km  (Figure  A7)  and  4  km  (not 
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shown),  respectively,  from  higher  resolution  NHMASS  output  of  potential  temperature 
and  vertical  velocity  are  consistent  with  previous  studies  of  inertia-gravity  wave  activity 
(Uccellini  and  Koch  1987).  The  wave  of  depression  at  16/2000  in  Figure  A6(a)  is  noted 
to  have  moved  to  the  Kentucky/West  Virginia  border  by  16/2200  [Figure  A6(b)],  along 
with  a  band  of  precipitation  out  of  phase  with  the  wave  (not  shown).  This  out-of-phase 
relationship  is  noted  to  be  one  piece  of  evidence  supporting  the  assessment  of  the  wave 
activity  as  IGW  activity  (Koch  and  Dorian  1988). 

A  second  test  for  IGW  activity,  quadrature  between  mid-tropospheric  vertical 
velocity  (w)  and  potential  temperature  (0)  fields,  is  evident  in  Figure  A8,  a  timeseries  of 
400-hPa  w  and  0  fields  at  a  point  in  northern  West  Virginia  from  the  FULL222 
simulation. 


Figure  A5.  NOAA  Global  Positioning  System  (GPS)  surface  pressure  trace  from 
Louisville  KY  (LOU)  GPS  site.  Valid  0000  UTC,  2/15/03  to  0000  UTC,  2/18/03.  Period 
of  IGW  activity  between  2000  UTC,  2/16/03  and  2300  UTC,  2/16/03  noted  with  black 
circle. 
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Figure  A6.  NHMASS  FULL18  (a)-(b)-(c)  one-hour  sea-level  pressure  tendency  (hPa/hr) 
for  2000  UTC,  2/16/03, 2100  UTC,  2/16/03,  and  2200  UTC,  2/16/03,  respectively.  Boxes 
in  (c)  denote  areas  over  which  covariances  of  perturbation  pressure  and  front  normal 
surface  wind  were  calculated  (see  text  for  details).  Hollow  circle  in  (a)  is  point  A  for 
timeseries  in  Figure  A8;  filled  circle  for  profiles  in  Figure  A9. 
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Figure  A7.  2-km  NHMASS  400-hPa  potential  temperature,  valid  0200  UTC,  2/17/03. 
Select  waves  noted  by  dashed  black  lines. 
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Figure  A8.  NHMASS  FULL2  time  series  of  (a)  400-hPa  vertical  velocity  [cm/s]  and  (b) 
400-hPa  potential  temperature  [K]  from  1830  UTC,  2/16/03  to  0330  UTC,  2/17/03  for 
point  A  displayed  in  Figure  A6(a). 


While  no  filtering  has  been  performed  on  the  analysis,  resulting  in  the  use  of  0 
rather  than  0’  and  w  rather  than  w’,  quadrature  between  400-hPa  w  and  0  is  clearly 
evident,  lending  credibility  to  the  assertion  of  the  wave  activity  as  IGW  activity.  An 
analysis  of  the  phase  relationship  between  surface  pressure  and  wave-normal  horizontal 
wind  speed  was  performed  (correlations  approaching  0.8)  and,  in  combination  with  the 
previous  tests  discussed,  supports  the  claim  of  inertia-gravity  waves  as  the  source  of  the 
mesoscale  banding  noted  in  the  available  datasets.  Due  to  the  limited  nature  of  asynoptic 
data  in  the  region  of  wave  propagation  and  lack  of  filtering  performed  in  the  above 
analyses,  no  definitive  claims  can  be  made.  As  a  result,  the  wave  activity  will  be  referred 
to  as  apparent  IGW  activity  within  this  paper. 
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A  TRANSITION  OF  IGWs  TO  TURBULENCE 

In  this  section,  the  process  through  which  vertically  propagating  apparent  IGW 
waves  underwent  steepening,  and  consequently  led  to  development  of  instability  and 
generation  of  turbulence  will  be  considered.  It  is  well  known  that  wave  breaking  is  often 
noted  in  regions  of  the  atmosphere  characterized  by  a  critical  layer  (where  wave  phase 
speed  is  equal  in  magnitude  to  the  mean  flow  in  the  direction  of  wave  propagation),  low 
Richardson  number,  and  a  wind  reversal  layer  (Jiang  and  Doyle  2004).  Near  the  critical 
level,  the  vertical  wavelength  of  the  waves  approach  zero  and  the  group  velocity  is 
directed  entirely  in  the  horizontal,  i.e.,  vertical  propagation  of  wave  energy  becomes  zero 
at  the  critical  level  (Lin  TBD).  The  orientation  of  gravity  waves  entirely  in  the  horizontal 
produces  fine-scale  instabilities  that  promote  wave  breaking  and  breakdown  of  the 
mesoscale  waves  into  turbulence. 

Considering  the  common  atmospheric  structure  in  the  vicinity  of  wave-breaking 
regions,  vertical  profiles  of  wave-front  normal  wind  speed  and  Richardson  number  from 
FULL222  simulations  will  be  considered.  Figure  A9(a),  a  profile  of  wave-front  normal 
wind  speed  as  determined  from  wave-front  positions  noted  in  Figure  A6,  indicates  a  deep 
layer  of  negative  wind  shear  above  approximately  180  hPa  with  multiple  layers  of 
enhanced  wind  shear.  Figure  A9(b)  indicates  that  a  low  Richardson  number  (less  than  5) 
is  simulated  between  about  120  hPa  and  80  hPa,  within  the  layer  of  negative  wind  shear. 
Although  Richardson  numbers  less  than  unity  have  been  observed  in  analyzing  critical 
layers  (Jiang  and  Doyle  2004;  Koch  and  Dorian  1988),  the  layer  of  interest  in  this  case  is 
in  the  very  stable  lower  stratosphere  (large  dO/dz),  where  the  Ri  calculation  is  likely 
skewed  larger  by  the  strong  static  stability,  despite  the  presence  of  large  values  of  vertical 
wind  shear. 

A  potential  range  of  critical  levels  has  been  included  in  Figure  A9(a),  determined 
from  both  calculations  of  phase  speed  from  linear  theory  and  a  subjective  analysis  of 
timeseries  of  70-hPa  velocity  divergence  fields  from  the  FULL222  simulation  (not 
shown).  Evident  in  a  comparison  of  Figures  A9(a)  and  A9(b)  is  the  fact  that  the  range  of 
possible  critical  levels  is  not  collocated  with  the  layer  of  the  low  Richardson  number.  A 
critical  limitation  of  this  analysis  is  that  a  total  wind  speed  (U+u’)  profile  has  been 
utilized,  rather  than  a  profile  of  background  wind  speed  (U).  It  has  been  noted  by  Fritts 
(1984)  that  a  critical  level  may  be  induced,  due  to  non-linear  effects  where  U+/-u’=c. 
While  this  raises  significant  questions  as  to  the  actual  location  of  a  critical  level  (or 
levels),  the  likelihood  of  turbulence  indicated  by  NHMASS  simulations  (to  be  discussed) 
and  the  presence  of  vertically  propagating  wave  activity  supports  the  presence  of  a  critical 
level  (or  levels)  within  the  lower  stratosphere. 
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Figure  A9.  0154  UTC,  2/17/03  222-m  NHMASS  profiles  of  (a)  wave-front  normal  wind 
speed  and  vectors  (m/s).  Range  of  critical  level  wind  speeds  (5-20  m/s)  noted  (see  text), 
(b)  Richardson  number,  vertical  profile  of  wave  front  normal  wind  (m/s)  included,  and  (c) 
Scorer  parameter  (solid  profile,  scaled  by  10^^  and  turbulence  kinetic  energy  dissipation 
(dashed  profile,  scaled  by  10^),  defined  in  text.  See  Figure  A6  for  location  of  profile. 
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The  vertical  variation  of  Scorer  parameter  (/^ ),  where 


£/"  U 


(A2) 


may  be  used  as  a  measure  of  the  potential  for  critical  level  wave  reflecting,  as  both  the 
Richardson  number  (Ri)  and  wind  reversal  requirements  discussed  above  are  accounted 
for  here  (Chun  1991).  This  is  true  since  Ri  is  proportional  to  static  stability,  which  can  be 
represented  by  N^,  the  Brunt  Vaisalla  frequency,  and  the  wind  reversal  layer  is  indicated 
by  the  curvature  of  u,  or  .  It  can  be  shown  from  linear  wave  theory  that  for  f  <  0,  an 
evanescent  flow  regime  exists  wherein  vertical  propagation  of  wave  energy  is  zero  and  for 
/^  >  0,  a  regime  exists  where  waves  are  free  to  propagate  vertically  (Lin  TBD). 

Therefore,  one  can  view  layers  where  [d(l^ )ldz  <  0  ]  as  layers  where  wave  modes 
transition  from  vertically  propagating  to  modes  where  the  vertical  wavelength  equals  zero. 
Such  a  transition  occurs  as  waves  approach  a  critical  level,  and  as  such  one  should  expect 
wave  steepening  and  potential  breaking/turbulence  generation  in  regions  where 
\d(V )jdz  <  0  ].  Thus,  from  the  Scorer  parameter  profile  in  Figure  A9(c),  the  layers  near 
100  hPa,  85  hPa,  50  hPa,  and  40  hPa  will  be  considered  for  wave  breaking  and  potential 
turbulence  generation. 

Steepening  wave  activity  can  be  seen  in  a  cross-sectional  analysis  of  vertical 
velocity  from  the  FULL222  simulation  at  17/0154  (Figure  AlO).  It  is  clear  that  the 
vertical  wavelength  of  the  waves  are  decreasing  with  increasing  height  in  the  lower 
stratosphere,  consistent  with  the  presence  of  a  critical  level  (or  levels).  One  can  see  that 
the  strongest  wave  is  located  near  150  hPa  in  the  center  of  the  cross-sectional  axis  with  an 
amplitude  of  greater  than  -0.9  m/s,  a  vertical  velocity  respectable  for  such  a  height  in  the 
atmosphere.  As  the  waves  steepen  and  largely  become  oriented  in  the  horizontal,  one 
expects  the  vertical  variation  of  potential  temperature  (0)  and  horizontal  wind  speed  (U)  to 
increase  as  the  alternating  phases  of  0  and  U  in  the  horizontal  are  tilted  into  the  vertical. 
These  changes  are  consistent  with  the  generation  of  mesoscale  instabilities,  two  of  which, 
convective  instability  and  Kelvin-Helmholtz  (or  shearing)  instability,  will  be  discussed 
next. 

Convective  instability,  present  where  dOejdz  <  0  (or  ddidz  <  0  for  very  dry  flow 
in  the  lower  stratosphere),  is  clearly  not  resolved  at  the  222-m  scale,  as  indicated  by  the 
lack  of  overturning  of  the  isentropes  in  Figure  All.  One  notes  dOfdz  >  0  everywhere 
within  the  lower  stratosphere.  Presumably,  with  finer  horizontal  and  vertical  grid 
resolution,  as  well  as  stronger  critical  levels  (i.e.,  more  wave  amplification),  overturning 
of  the  isentropes  in  the  lower  stratosphere  would  be  explicitly  resolved.  A  second 
instability,  Kelvin-Helmholtz  instability,  produces  growing  disturbances  through 
extraction  of  kinetic  energy  present  in  strongly  sheared  flow  (|5f//5z| »  0 )  for  use  in 

raising  or  suppressing  stable  fluid  (A^^  >0).  This  instability  is  most  effective  for  strong 
shear  and  weak  static  stability  (i.e.,  Ri  <  14)  (Lin  TBD). 
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While  multiple  layers  of  enhanced  vertical  wind  shear  are  apparent  in  a  lower- 
stratospheric  cross-section  of  total  wind  speed  (Figure  Al  1),  the  lack  of  layers  with  both 
strong  enough  shear  and  sufficiently  weak  static  stability  precludes  Kelvin-Helmholtz 
instabilities  from  developing  in  the  FULL222  simulation.  It  must  be  emphasized  here  that 
the  convective  and  orographic  forcing  supporting  the  wave  breaking  layers  was  weak 
during  this  case,  increasing  the  difficulty  of  simulating  full  wave  breaking.  NHMASS 
simulations  by  another  NHMASS  modeler  (Paul  Suffern)  of  a  separate  case  study  of 
inertia-gravity  wave  breaking  in  a  region  of  stronger  convective  forcing  with  222-m  grid 
spacing  have  explicitly  resolved  the  above  instabilities  and  through  use  of  a  71 -m 
resolution  simulation,  turbulent  mixing,  indicating  the  mesoscale  model’s  ability  to 
resolve  such  features.  Lane  et  al.  (2003),  in  a  numerical  modeling  study  of  turbulence 
above  deep  convection,  was  able  to  resolve  regions  of  turbulence  with  a  3-dimensional 
cloud  model  with  grid  spacing  approximately  'A  that  used  in  the  finest  of  the  simulations 
in  this  study.  The  smallest  resolvable  feature  in  this  simulation,  on  the  order  of  888  m 
(4Ax  wave),  is  much  larger  than  the  scale  of  most  turbulent  eddies  (~100-  to  300-m). 


Figure  A 10.  222-m  NHMASS  vertical  cross-section  of  vertical  velocity  (cm/s),  solid 
(positive),  dashed  (negative),  valid  0154  UTC,  2/17/03.  See  inset  for  cross-section  axis. 
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Figure  All.  222-m  NHMASS  vertical  cross-section  of  total  wind  speed  (m/s)  (shaded, 
solid  contour)  and  potential  temperature  (K)  (dashed),  valid  17/0154.  Regions  of  potential 
Kelvin-Helmholtz  instability  indicated  by  ovals.  For  cross-section  axis  see  inset  in  Figure 
AlO. 


Caution  must  be  exercised  in  interpreting  mesoscale  model  output  at  such  fine 
horizontal  resolutions  (such  as  71  m),  as  the  unequal  vertical  and  horizontal  resolutions 
negates  the  proper  simulation  of  large  turbulent  eddies.  Additionally,  due  to  the  need  for 
increased  vertical  resolution  both  in  the  boundary  layer  and  in  the  lower  stratosphere,  a 
discontinuity  in  vertical  resolution  exists  in  the  lower  stratosphere  near  110  hPa.  This 
discontinuity,  on  the  order  of  a  change  of  vertical  resolution  of  ~500  m/10  hPa,  raises 
doubt  as  to  the  degree  of  simulated  wave-breaking,  as  some  wave  reflection  in  the 
NHMASS  simulations  near  1 10  hPa  may  be  numerical  in  origin.  The  use  of  finer  scale 
simulations  in  this  study  was  prevented  by  extreme  near-surface  instabilities  generated 
due  to  errors  in  the  boundary  layer  scheme  in  the  vicinity  of  fine-scale  terrain. 

Despite  the  insufficient  resolution  with  which  to  explicitly  resolve  the  above 
instabilities  (and  therefore  turbulent  mixing),  the  resolved  fields  within  the  FULL222 
simulation  do  indicate  the  likelihood  of  sub-grid  scale  turbulence.  The  first  large-scale 
measure  of  turbulence  to  be  considered,  the  NCSUl  index,  defined  as 


Ncsm^^  •VU 


(A3) 
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has  been  shown  to  be  sensitive  to  regions  of  strong  anticylonic  shear  in  the  presence  of 
low  Ri  (conditions  conducive  for  clear-air  turbulence  [CAT])  (Kaplan  et  al.  2004). 
Inertia-gravity  wave  activity  perturbs  the  large-scale  flow,  producing  the  regions  of  strong 
lateral  shear  and  weak  dynamic  stability  that  the  NCSUl  parameter  is  designed  to  be 
sensitive  to.  Figure  A 12  indicates  a  region  of  enhanced  NCSUl  values  moving  northeast 
toward  the  location  of  the  aforementioned  vertical  profiles  and  cross-sections  that  were 
suggestive  of  wave  breaking.  Of  interest  is  the  fact  that  the  region  of  high  NCSUl  index 
values,  apparently  maximized  at  17/0154  [Figure  A 1 2(d)]  is  positioned  immediately  west 
of  a  roughly  150-m  steep  terrain  ridge  across  northeastern  West  Virginia  (Figure  A13). 
One  possibility  is  that  mountain-wave  drag  due  to  low-level  easterly  flow  normal  to  the 
ridge  induces  a  critical  level  (or  levels)  in  the  tower  stratosphere  (Lin  TBD),  The  apparent 
inertia-gravity  waves,  propagating  toward  northeastern  West  Virginia,  encounter  these 
mountain- wave  induced  critical  levels,  steepen,  and  break-down  into  turbulence. 


Figure  A12.  222-m  NHMASS  NCSUl  index  [scaled  by  lO’^]  valid  (a)  0145  UTC, 
2/17/03,  (b)  0148  UTC,  2/17/03,  (c)  0151  UTC,  2/17/03,  and  (d)  0154  UTC,  2/17/03.  See 
text  for  definition  of  NCSUl  index. 
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Figure  A13.  222 -m  NHMASS  model  terrain  [m]. 


The  mean  turbulence  kinetic  energy  (TKE),  defined  as 


TKE 

m 


+v'^ 


(A4) 


(Stull  1988),  is  the  most  direct  measure  of  turbulence  available,  A  turbulence  budget 
equation  can  be  written  for  mean  turbulence  kinetic  energy, 
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assuming  horizontal  homogeneity  (Stull  1988).  In  Equation  (A5),  term  I  is  the  TKE 
tendency,  term  II  and  III  are  shear  production/loss  terms,  term  IV  is  the  buoyant 
production/dissipation  term,  term  V  represents  vertical  diffusion,  and  term  VI  represents 
the  viscous  dissipation  of  TKE.  The  boundary  layer  parameterization  used  for  all 
NHMASS  simulations  in  this  study  is  based  on  the  1.5  order  TKE  closure  scheme 
developed  by  Therry  and  Lacarrerre  (1983),  wherein  a  form  of  the  TKE  budget  Equation 
(A5)  is  used  in  closing  the  model  equations.  Post-processing  of  NHMASS  output  makes 
use  of  the  TKE  scheme  in  generating  TKE  profiles,  yet  due  to  the  apparently  weak  nature 
of  the  turbulence  in  the  lower  stratosphere,  the  output  of  raw  TKE  proved  to  be  of  little 
value  as  small  quantities  of  TKE  are  output  as  zero.  Rather,  TKE  dissipation  (term  VI  in 
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Equation  (A5))  calculations  were  relied  upon  in  this  study  and  are  generated  using  the 
empirical  relationship, 

(A6) 

LK 

where  LE  is  the  dissipation  length  (Kaplan  et  al.  2000).  While  a  more  rigorous  analysis 
would  involve  calculating  terms  I-V  in  Equation  (A5)  and  estimating  f  as  a  remainder, 
this  empirical  method  does  appear  to  have  some  value. 

Considering  a  profile  of  s  [dashed  profile  in  Figure  A9(c)],  one  notes  a  good 
correlation  between  the  layers  of  possible  critical  level-induced  wave  breaking  (1^ 
decreasing  with  height)  and  the  layers  of  enhanced  e .  One  also  notices  that  the  greatest 
values  of  e  are  not  associated  with  the  more  pronounced  layers  of  5/5z(/^)<  0  in  the  80- 

100-hPa  layer  but,  rather,  are  maximized  near  the  50-hPa  level  where  gradients  of  P  are 
smaller.  One  possibility  for  this  verticals^  structure  is  that  with  each  successive  critical 
level  encountered  by  vertically  propagating  waves,  each  level  is  not  strong  enough  to  limit 
further  vertical  movement  of  wave  energy,  and  so  multiple  wave-breaking  episodes  occur 
with  the  ^SO-hPa  level  receiving  the  benefits  of  the  subsequent  turbulence  generation. 
Nevertheless,  the  two  vertical  profiles  [Figure  A9(c)]  and  the  NCSUl  analysis  (Figure 
A 12)  appear  to  support  the  idea  that  (1)  vertical  variation  in  Scorer  parameter  can  indicate 
possible  wave-breaking  layers  and  (2)  the  variations  in  the  mesoscale  horizontal  fields 
represented  by  the  enhanced  values  of  the  NCSU 1  index  at  50  hPa  were  a  reasonable 
indicator  of  50-hPa  turbulence  potential  across  northeastern  West  Virginia. 

5.  SUMMARY  AND  CONCLUSIONS 

The  threat  of  atmospheric  turbulence  to  the  aviation  community  motivated  this 
study  of  inertia-gravity  wave  breaking  and  generation  of  lower-stratospheric  turbulence 
during  the  2003  Presidents’  Day  winter  storm.  The  ability  of  a  mesoscale  model,  the 
Non-hydrostatic  Mesoscale  Atmospheric  Simulation  System  (NHMASS)  version  6.3,  to 
simulate  the  development  of  apparent  inertia-gravity  wave  activity  in  the  mid-troposphere 
and  subsequent  transition  to  turbulence  in  the  lower  stratosphere  has  been  evaluated. 
Simulation  of  both  the  pre-existing  unbalanced  mid-tropospheric  flow  and  the  resultant 
apparent  inertia-gravity  wave  activity  have  been  shown  to  be  consistent  with  the  synoptic 
and  limited  asynoptic  data  available.  The  vertically  propagating  wave  activity  was  noted 
to  propagate  toward  northeastern  West  Virginia,  losing  coherency  north  and  east  of  the 
region  (not  shown).  The  Scorer  parameter,  a  measure  of  both  vertical  decrease  of  static 
stability  and  the  presence  of  negative  vertical  wind  shear  (both  consistent  with  critical 
layer  theory),  was  evaluated  over  northeastern  West  Virginia.  This  analysis  of  the  vertical 
structure  of  the  upper  troposphere  and  lower  stratosphere  revealed  a  number  of  potential 
levels  where  vertical  propagation  of  wave  energy  was  limited  and  the  likelihood  of  wave 
breaking  increased. 

Cross-sections  of  vertical  velocity  indicated  vertically  propagating  wave  activity  in 
the  upper  troposphere,  with  vertical  wavelengths  decreasing  with  height  in  the  lower 
stratosphere,  consistent  with  the  presence  of  a  critical  level  (or  levels)  in  that  layer.  Two 
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mesoscale  instabilities  noted  to  develop  in  regions  of  breaking  inertia-gravity  wave 
activity,  convective  instability  and  Kelvin-Helmholtz  (or  shearing)  instability  (Lane  et  al 
2003),  were  not  resolved  with  222-m  horizontal  resolution,  yet  the  existence  of  sub-grid 
scale  turbulence  was  inferred  from  two  turbulence  indicators,  the  NCSU 1  index  and 
turbulence  kinetic  energy  (TKE)  dissipation.  The  NCSUl  index,  sensitive  to  horizontal 
gradients  of  momentum  and  small  values  of  the  Richardson  number,  was  evaluated  at  50 
hPa  (one  of  the  layers  indicated  by  Scorer  parameter  profiles  as  a  possible  wave-breaking 
layer),  revealed  the  likelihood  of  turbulence  in  northeastern  West  Virginia.  Bands  of  high 
NCSUl  index  were  noted  to  have  been  positioned  along,  and  in  the  lee  of,  100-  to  300-m 
tall  terrain  ridges  in  the  region,  where  mountain  wave  activity  potentially  induced  the 
critical  levels  deduced  in  the  aforementioned  vertical  profiles.  A  comparison  of  vertical 
profiles  of  TKE  dissipation  and  Scorer  parameter  indicated  a  high  correlation  between  the 
potential  wave-breaking  layers  indicated  by  the  vertical  variation  of  Scorer  parameter  and 
the  regions  of  enhanced  TKE  dissipation. 

This  study  indicates  the  suitability  of  two  parameters  in  detecting  the  potential  for 
wave  breaking  and  lower-stratospheric  turbulence  utilizing  simulations  with  insufficient 
resolution  with  which  to  explicitly  resolve  either  full  wave  breaking  or  turbulent  mixing. 
One  parameter  is  sensitive  to  lateral  changes  associated  with  the  increased  anticyclonic 
shear  produced  by  amplifying  inertia-gravity  waves  (NCSUl),  and  the  second  is  sensitive 
to  vertical  variation  of  the  perturbed  flow  suggestive  of  wave-breaking  potential  (Scorer 
parameter).  Improving  the  predictability  of  turbulence  is  an  effort  sure  to  have  many 
positive  effects  on  the  aviation  community  and  the  public-at-large. 

6.  FUTURE  WORK 

Additional  case  studies  are  currently  being  analyzed,  applying  the  same 
methodology  described  here,  in  order  to  further  test  the  predictability  of  turbulence  in  the 
lower  stratosphere  using  the  NCSU  1  index  and  vertical  variation  of  Scorer  parameter. 
Lastly,  the  possibility  of  utilizing  high  temporal  and  spatial  resolution  NOAA-16 
Advanced  Microwave  Sounding  Unit  (AMSU)  brightness  temperature  data  to  better 
validate  the  existence  of  inertia-gravity  wave  activity  in  NHMASS  simulations  is  being 
pursued. 
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Appendix  B 

A  TKE  Budget  Analysis  in  the  Upper  Troposphere/Lower  Stratosphere  for  a 
Vertically  Propagating  Gravity  Wave 


1.  INTRODUCTION 

The  application  of  turbulent  kinetic  energy  (TKE)  and  the  TKE  tendency  equation 
are  usually  applied  to  quantify  turbulence  and  turbulent  tendency  within  the  planetary 
boundary  layer  (PBL).  The  importance  of  TKE  and  TKE  tendency  is  far  greater  within 
the  PBL  than  the  free  atmosphere  on  average.  However,  turbulence  is  of  great  importance 
in  the  upper  atmosphere,  especially  to  aviation  interests. 

The  United  States  Air  Force  has  commissioned  a  project  with  the  North  Carolina 
State  University’s  Mesoscale  Modeling  and  Dynamics  Lab  to  investigate  the  role  of 
convection  on  generating  turbulence  in  the  lower  stratosphere.  Several  different  case 
studies,  representing  various  different  synoptic  and  mesoscale  environments  and  different 
convective  modes,  are  to  be  simulated  using  numerical  models.  The  importance  of 
vertically  propagating  gravity  waves  near  the  tropopause  would  have  great  influence  on 
the  turbulence  generated  in  these  layers.  Aircraft  who  choose  to  fly  over  convection, 
especially  Air  Force  spy  planes  that  fly  at  high  levels,  may  unexpectedly  experience 
turbulence  at  vertical  levels  far  above  the  convection. 

A  TKE  budget  analysis,  involving  the  TKE  tendency  equation,  is  done  explicitly 
for  the  12  December  2002  large  amplitude  gravity  wave  over  far  eastern  Texas  and 
western  Louisiana.  Figure  B 1  shows  the  stratospheric  wave  for  the  model  simulation 
described  in  this  paper  as  a  cross-section  across  the  region.  Though  rudimentary,  this  first 
attempt  at  a  TKE  budget  analysis  is  aimed  at  providing  an  initial  insight  and  motivation 
for  future  research  related  to  the  author’s  thesis  and  the  Air  Force  project’s  other  case 
study  simulations.  The  results  presented  within  this  paper  in  no  way  represent  a  truly 
accurate  and  rigorous  TKE  budget  analysis,  but  provide  a  good  qualitative  starting  point 
for  future  research  related  to  TKE  and  eddy  dissipation  (EDR)  budget  calculations. 
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Figure  Bl.  Cross-section  showing  the  total  wind  field  (color)  and  the  isentropes 
(contoured)  for  the  222-m  horizontal  STRATONHMASS  simulation  of  the  12  December 
2002  large  amplitude  gravity  wave  valid  at  2227  UTC,  12  December. 


2.  METHODOLOGY 
2.1.  Reynolds  Decomposition 

For  the  purposes  of  the  TKE  tendency  equation,  it  is  useful  to  employ  Reynolds’ 
decomposition.  For  any  scalar  quantity,  its  value  at  any  time  step  or  position  can  be 
expressed  as  a  combination  of  an  average/  basic  state  and  a  perturbation  from  that  basic 
state,  namely: 


n  =  Q  +  0'  (Bl) 

For  the  output  from  a  numerical  model,  the  arbitrary  value  Q  would  be  the  output  at  each 
interval  specified  by  the  user.  The  calculation  of  this  arbitrary  variable  H  is  inherently  an 
average  itself,  coming  from  the  model’s  algorithm  using  the  values  at  each  time  step.  For 
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this  study,  the  parameter  Q  is  considered  to  be  representative  of  a  true  time  step  value,  and 
not  an  average.  The  complications  that  arise  from  needing  to  modify  the  computational 
scheme  to  extract  the  value  of  Q  at  each  time  step  are  beyond  the  current  study. 

If  is  known  and  is  treated  as  a  set  of  model  simulated  values  over  a  given  time 

range,  the  value  of  the  basic  state/average  Q  can  be  calculated  by  simply  taking  a  time 
average,  such  that: 


—  (f2.  +  Ci^  ^3  “t"  ...fij) 

Q  - L_i - '1  (B2) 

i 

where  i  is  the  number  of  output  times  the  variable  Q  is  output.  The  user  can  choose 
different  ranges  to  calculate  the  average  field  over.  Because  of  time  restraints,  an  average 
over  the  entire  model  run  is  done  for  the  case  study  presented  herein.  If  O  is  known,  then 
the  perturbation  can  be  calculated  by  solving  for  the  perturbation  value  in  Equation  (Bl): 

Q'  =  0-n  (B3) 

Though  this  method  is  rudimentary,  it  is  most  convenient  for  the  data  types  presented  in 
the  Air  Force  case  studies  and  provide  a  good  qualitative  representation  of  the  average  and 
perturbation  fields  for  any  quantity  from  the  model  output.  As  aforementioned,  the  user 
selects  the  amount  of  time  the  variable  Q  is  averaged  and,  thus,  selects  the  time  scale  at 
which  the  perturbations  are  evaluated.  To  capture  small-scale  features,  it  is  advantageous 
to  select  a  shorter  time  range  to  calculate  the  basic  state  and,  henceforth,  capture  the 
smaller  scale  perturbations  in  the  Cl  field.  If  the  length  of  the  actual  model  run  is  short, 
and  the  horizontal  and  vertical  grid  spacing  are  small  enough  to  capture  the  features  of 
interest,  this  restriction  can  be  relaxed. 

2.2.  TKE  Tendency  Equation 

In  order  to  quantify  and  understand  the  contributions  to  turbulence  production,  an 
analysis  of  each  term  in  the  TKE  tendency  equation  is  needed.  Numerical  models,  such  as 
the  Non-Hydrostatic  Stratospheric  Mesoscale  Atmospheric  Simulation  System 
(STRATONHMASS)  output  TKE  at  each  specified  output  time.  The  TKE,  defined  as: 


e  +w'^  (B4) 

is  calculated  at  each  time  step  and  averaged  for  each  output  time.  The  model  user  is 
unable  to  see  variances  change  from  each  model  time  step.  In  addition,  the  TKE  scheme 
does  not  explicitly  show  or  calculate  the  various  terms  in  the  TKE  tendency  equation  for 
mesoscale  numerical  models. 
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From  Stull  (1988),  the  TKE  tendency  equation  can  be  written  as: 


de  -Be  -Be  —Be  g  ^  , 

—  =  ~u - V - w — +^(w  0^  ) 

Bt  Bx  By  Bz  9^ 


1  B 


—  .Bu  B 


oz  oz  p  oz 


(B5) 


Equation  (B4)  is  valid  if  the  model  coordinate  system  is  aligned  with  the  mean  wind,  but 
does  not  necessitate  the  need  for  horizontal  homogeneity  as  prescribed  by  Stull  (1988, 
hereafter  S88).  The  only  term  on  the  left-hand  side  of  Equation  (B5)  is  denoted  the 
storage  term  by  S88  and,  for  the  duration  of  this  paper,  is  referred  to  as  the  TKE  tendency. 
The  TKE  tendency  can  be  computed  over  a  given  time  range  if  the  variables  on  the  right- 
hand  side  were  known  over  the  same  time  range. 

The  first  term  on  the  right-hand  side  of  Equation  (B5)  is  the  creation/destruction  of 
turbulence  by  the  advection  of  the  mean-u  wind,  which  will  be  denoted  Term  1:  U- 
Advection.  Similarly,  the  second  and  third  terms  are  the  creation/destruction  of  turbulence 
by  the  advection  from  the  mean-v  and  mean-w  wind,  and  are  denoted  Term  2:  V- 
Advection  and  Term  3:  W-Advection.  Authors  such  as  S88  neglect  the  horizontal 
advection  by  assuming  the  horizontal  flow  is  homogenous,  and  are  much  smaller  in 
comparison  to  Term  3:  W-Advection  for  situations  where  turbulence  is  important.  Due  to 
the  rapid  fluctuations  in  u  and  v  during  the  case  study  presented  in  Section  3,  it  was 
necessary  to  investigate  the  overall  contributions  from  Term  1  and  Term  2  to  the  overall 
TKE  tendency  profile.  Term  3,  which  includes  the  mean  vertical  motion  field,  is 
hypothesized  to  be  important  in  locations  where  large  variations  in  TKE  are  overlapped 
with  strong  background  vertical  motion  from  convective  motions. 

The  fourth  term  on  the  right-hand  side  of  Equation  (B5)  is  referred  to  as  the 
buoyancy  production/consumption  term  in  S88,  and  is  referred  to  hereafter  as  Term  4: 
Buoyancy.  The  physical  interpretation  of  this  term  is  that  positively  buoyant  (negatively 
buoyant)  parcels  act  to  produce  (consume)  TKE.  The  product  of  the  perturbations,  which 
are  the  vertical  velocity  and  the  virtual  potential  temperature,  mimic  a  vertical  heat  flux 
term.  S88  states  Term  4’s  biggest  contribution  within  the  PBL  comes  on  days  where  free 
convection  and  subsequent  strong  surface  heating  takes  place.  For  stratospheric 
applications,  it  is  hypothesized  that  areas  where  perturbations  in  vertical  motion  and 
virtual  potential  temperature  from  convection  at  and  below  the  tropopause  would  produce 
large  and  small  variations  in  Term  4. 

The  fifth  term  on  the  right-hand  side  of  Equation  (B5)  is  known  as  the  mechanical 
shear  production  term,  from  S88,  and  is  referred  to  hereafter  as  Term  5:  Mechanical 
Shear.  According  to  S88,  when  the  turbulent  momentum  flux  is  in  the  presence  of  a  large 
mean  wind  shear,  the  two  interact  positively  to  generate  more  turbulence.  The  product  of 
the  perturbations  of  the  u  and  w  wind  is  a  momentum  flux,  and  this  flux  is  multiplied  by 
the  basic  state  mean  wind  shear.  Term  5  would,  in  theory,  be  maximized  for  stratospheric 
applications  where  a  large  mean  wind  shear  is  found  over  the  averaging  time  period.  As 
will  be  demonstrated  for  the  12  December  2002  case,  relatively  small  vertical  propagating 


46 


waves  are  capable  of  producing  large  mean  wind  shear  values  around  their  periphery.  The 
large  wind  shear  over  a  short  distance  would  act  to  produce  a  non-negligible  contribution 
to  the  TKE  tendency. 

The  sixth  term  of  Equation  (B5)  represents  the  contribution  of  TKE  from  the  rate 
of  change  in  the  vertical  of  turbulent  transport,  hereafter  referred  to  as  Term  6:  Turbulent 
Transport.  The  term  that  results  is  similar  to  a  vertical  TKE  flux;  and  the  contribution  to 
TKE  tendency  comes  from  the  vertical  derivative.  According  to  S88,  on  a  local  scale, 
Term  6  acts  to  either  enhance  (retard)  TKE  generation  depending  on  where  there  is  flux 
convergence  (divergence).  If  there  is  more  flux  into  a  given  layer  than  leaves  in  that 
layer,  then  the  magnitude  of  the  TKE  increases  (S88).  Stull  (1988)  comments  that  Term  6 
really  does  not  create  or  destroy  TKE,  but  rather  re-distributes  it  within  a  given  layer.  For 
a  stratospheric  case,  where  vertical  perturbations  are  large,  the  author  hypothesizes  that 
the  greatest  contribution  to  TKE  creation  (destruction)  would  be  in  a  situation  where  large 
values  of  TKE  are  juxtaposed  with  strong  upward  (downward)  perturbation  velocities. 

The  seventh  term  of  Equation  (B5)  represents  the  creation  and  destruction  of  TKE 
by  pressure  perturbations,  which  is  called  the  “pressure  correlation”  term  by  S88. 
Hereafter,  the  pressure  correlation  term  is  referred  to  as  Term  1:  Pressure  Perturbation. 
S88  comments  that  very  little  is  known  about  this  term  from  observational  data  since  it  is 
difficult  to  measure  a  static  pressure  perturbation  flux  in  the  atmosphere.  However,  in  a 
numerical  model  such  as  the  STRATONHMASS,  the  pressure  perturbation  term  can  be 
output  directly  from  the  model.  Term  7  can  then  be  calculated  directly,  given  that  the 
TKE  is  output  at  the  same  time  interval.  The  author  hypothesizes  that  the  contribution 
from  Term  7  would  be  largest  where  pressure  perturbations  and  TKE  are  juxtaposed,  most 
likely  above  strong  convective  updrafts,  where  TKE  is  large  due  to  momentum  variances. 

The  final  term  in  Equation  (B5)  is  the  TKE  dissipation,  which  is  always  negative 
and  acts  to  reduce  TKE.  The  eighth  term  is  referred  hereafter  as  Term  8:  Dissipation. 
Term  8  is  also  referred  to  as  the  eddy  dissipation  rate  (EDR).  The  dissipation  is  from 
molecular  destruction  and,  thus,  an  empirical  scheme  must  be  employed,  as  it  cannot  be 
calculated  directly  from  the  model  output  as  the  other  terms  in  Equation  (B5).  The 
STRATONHMASS  outputs  the  dissipation  at  each  output  along  with  the  rest  of  the 
parameters  listed  in  Equation  (B4).  The  TKE  dissipation  (EDR)  scheme  is  not  described 
here,  but  is  a  high  order  scheme  that  does  output  the  dissipation  at  every  level  outside  the 
PBL. 

2.3.  TKE  Tendency  Equation  in  Pressure  Coordinates 

Though  not  described  explicitly  in  S88,  for  a  model  data  set,  it  is  often  convenient 
to  express  variables  in  pressure  coordinates  as  opposed  to  height  coordinates.  For  most 
data  sets  from  field  studies,  parameters  are  measured  as  a  function  of  height,  which  is 
known  given  the  height  of  the  instrument.  However,  in  a  numerical  model,  the  vertical 
axis  is  most  often  expressed  as  a  normalized  coordinate  system,  such  as  a,  which  can 
easily  be  output  into  corresponding  pressure  coordinates.  Graphical  interpretation 
packages  such  as  GEMPAK  most  often  work  with  data  that  is  plotted  in  pressure 
coordinates. 
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In  addition  to  the  overall  general  use  of  pressure  versus  height  coordinates,  the 
height  spacing  between  vertical  a  levels  is  uneven  for  the  STRATONHMASS  model. 

The  model’s  post-processing  package  automatically  extrapolates  the  a  levels  into 
corresponding  pressure  values,  but  does  not  for  height  coordinates.  The  uneven  nature  of 
the  height  spacing  between  a  levels  introduces  a  major  data  analysis  problem  that  may 
require  quite  a  bit  of  programming  to  address.  Therefore,  for  the  purposes  of  this  study, 
the  use  of  pressure  coordinates  is  most  convenient. 

For  the  purposes  of  this  study  using  a  STRATONHMASS  simulation-produced 
data  set.  Equation  (B5)  must  be  expressed  in  pressure  coordinates.  To  do  so,  we  consider 
the  vertical  derivative  of  the  arbitrary  quantity  Q  expressed  as: 


50  dp 
dz  dp  dz 


(B6) 


from  a  generalized  coordinate  transformation,  taken  from  Bluestein  (1993).  If  we  assume 
that  the  atmosphere  is  hydrostatically  balanced  within  the  layer  of  the  derivative,  then  we 
can  invoke  the  hydrostatic  approximation,  such  that  Equation  (B6)  becomes: 


50  50 

dz  dp 


(B7) 


The  author  acknowledges  that  the  output  used  is  from  a  non-hydrostatic  model,  and  that 
the  hydrostatic  balance  is  not  applicable  to  convective  systems.  However,  the  hydrostatic 
approximation  is  only  made  over  the  difference  between  the  levels  in  which  the  derivative 
is  calculated  (see  Section  2.4).  The  pressure  differences  between  the  levels  are  between  5 
and  10  millibars,  and  even  for  the  stratosphere,  only  correlate  to  a  0.2-km  difference  in 
levels.  Thus,  the  hydrostatic  approximation  is  likely  valid  within  those  layers  for  the 
duration  of  the  model  run.  Sensitivity  studies  to  a  higher-order  coordinate  transformation 
would  be  needed  to  exactly  quantify  the  error  introduced  by  using  the  hydrostatic 
approximation  with  the  levels  between  the  derivatives. 

Using  Equation  (B7),  we  can  rewrite  Equation  (B5)  as: 


5e 

¥ 


-de  -de  —  de 

acc-u - V —  +  wpg  — 

dx  dy  dp 


+ 


0. 


+  ( u'w')pg^  +  pg~( w'e)  +  g-^(w'p’)  -  s 
dp  dp  dp 


(B8) 


Equation  (B8)  is  the  TKE  tendency  equation  written  in  pressure  coordinates,  and  founds 
the  basis  for  the  pressure-based  TKE  budget  presented  in  Section  3.  For  the  simplicity  of 
the  calculations,  the  author  takes  the  average  air  density  to  be  1.25  kg/m^  and  gravity  as 
9.80665  m/s^ 
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2.4.  Finite  Difference  Scheme 

In  order  to  represent  Equation  (B8)  to  calculate  derivatives  using  actual  data,  finite 
difference  schemes  must  be  implemented  to  approximate  the  derivatives.  Following  Lin 
(2005),  the  most  stable  and  useful  finite  difference  scheme  is  a  centered  in  time  and 
centered  in  space  approach.  From  Equation  (B8),  we  see  that  only  spatial  derivatives  are 
necessary  to  quantify  the  budget,  since  we  are  solving  for  a  time  tendency  on  the  right 
hand  side. 

The  horizontal  derivatives,  seen  in  Term  1  and  Term  2  of  Equation  (B8),  can  be 
calculated  directly  using  the  GEMPAK  numerical  graphics  package.  Due  to  time 
constraints,  the  horizontal  derivatives  for  the  case  study  presented  in  Section  3  are 
calculated  using  GEMPAK’s  horizontal  derivative  scheme.  The  functions  DDX  and  DDY 
are  valid  for  any  scalar  field  and  follow  a  second-order  centered  spatial  finite 
difference. 

Because  the  choice  of  vertical  axis  is  left  to  the  user,  the  GEMPAK  numerical 
graphics  package  does  not  calculate  a  vertical  derivative  of  any  sort.  From  Equation  (B8), 
we  see  that  Terms  3,  5, 6,  and  7  all  require  a  vertical  derivative  of  some  scalar  quantity  or 
product  of  quantities.  In  order  to  account  for  these  terms,  one  must  develop  a  finite 
difference  scheme  that  is  consistent  with  the  horizontal  derivatives  calculated  in  Term  1 
and  Term  2.  From  Lin  (2005),  a  first-order  centered  spatial  finite  difference  for  some 
arbitrary  quantity  Q  is  done  vertically,  such  as  shown  in  Equation  (B9): 


dp  2Af> 


(B9) 


The  scalar  quantity  Q  is  taken  at  time  step  t  from  two  different  pressure  levels  above  and 
below  the  level  at  which  the  derivative  is  calculated.  Since  pressure  has  the  opposite  sense 
as  height,  it  must  be  pointed  out  that  the  first  observation  corresponds  to  the  value  of  Q  at 
a  higher-pressure  level,  meaning  lower  height.  Fortran  programs  were  written  to  calculate 
the  vertical  derivaties  for  the  various  scalar  quantities  shown  in  Equation  (B8),  and  are 
described  in  further  detail  in  Section  3. 

Though  the  horizontal  derivatives  are  of  higher  order  accuracy  than  Equation  (B9), 
the  two  schemes  are  still  comparable  from  a  semi-quantitative  viewpoint.  Due  to  time 
constraints,  a  comparison  of  Equation  (B9)  to  a  second-order  centered  vertical  scheme 
was  unable  to  be  completed,  but  is  planned  to  be  part  of  future  work  on  the  project. 

3.  TKE  BUDGET  ANALYSIS:  12  DECEMBER  2002  STRATONHMASS 
SIMULATION  IN  PRESSURE  COORDINATES 

3.1.  Model  Simulation  Description 

The  STRATONHMASS  model  was  run  at  a  horizontal  grid  spacing  of  Ax,  Ay  = 
222  m  for  the  12  December  2002  case  using  the  pressure  coordinated  transformation 
discussed  in  Section  2.3.  The  model  was  initialized  at  2209  UTC,  12  December,  and  ran 
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approximately  35  minutes,  with  the  last  model  output  time  of  2244  UTC,  12  December. 
The  model  was  initialized  using  NCAR  NARR  reanalysis  data  that  includes  data  above 
100  millibars,  making  computations  in  the  stratosphere  possible. 

Figure  B2  shows  the  terrain  map  for  this  222-meter  simulation.  The  grid  area  is 
very  small  and  only  covers  a  small  portion  of  far  eastern  Texas  and  western  Louisiana. 

The  simulation  was  run  using  the  Kain-Fritch  convective  parameterization  scheme  and 
Lin  et  al.  (1983)  microphysics  package.  The  model  was  run  with  a  TKE  option  for  the 
PBL,  but  this  is  of  no  concern  since  we  are  dealing  with  the  upper  level  regions  for  the 
TKE  budget. 

The  model’s  time  step  is  approximately  0.07  s,  but  the  output  from  the  model  run 
was  only  available  every  54  seconds.  This  totaled  40  model  simulation  times  for  each 
field  over  the  35-minute  model  run.  Due  to  a  time  formatting  error  within  the  GEMPAK 
graphics  package,  five  of  these  40  model  simulations  times  had  to  be  neglected.  Thus,  the 
TKE  budget  is  performed  using  35  total  model  simulation  times  and  the  initial  state. 

The  model  was  run  with  90  vertical  ct  levels.  For  the  purposes  of  the  TKE  budget 
analysis,  however,  only  data  between  300  mb  and  10  mb  was  considered.  While  the  entire 
vertical  structure  of  the  model  can  be  ingested  into  the  programs  and  scripts  necessary  to 
run  the  TKE  budget,  the  computational  time  required  is  not  feasible  if  one  is  only 
interested  in  stratospheric  applications.  Model  data  that  lie  between  the  300  millibar  and 
100  millibar  surfaces  are  available  every  10  millibars,  while  data  above  100  millibars  to 
the  model  top  at  10  millibars  are  available  every  5  millibars. 
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Figure  B2.  Domain  terrain  map  for  the  222-m  STRATONHMASS  simulation,  shown  in 
meters,  for  the  12  December  2002  case. 


In  the  model,  as  was  true  with  satellite  and  ground-based  GPS  observations,  strong 
convection  forced  a  large-amplitude  gravity  wave  that  propagated  within  a  wave  duct 
north  of  a  strong  warm  front/stable  layer.  In  Paul  Suffem’s  222-m  output  data  available 
on  the  Mesolab  server,  a  strong  vertically  propagating  gravity  wave,  initiated  by  the 
convection  below,  traveled  up  the  troposphere  and  into  the  lower  stratosphere.  The 
convection  top/overshoot  from  the  w-relative  wind  fields  was  around  170  mb,  penetrating 
50  mb  into  the  stratosphere.  According  to  model-derived  soundings,  the  tropopause  was 
iocated  around  210  mb.  From  the  plots  of  the  total  wind  fields  and  isentropes  (not 
shown),  the  levels  at  which  the  budget  are  most  interesting  are  between  120  mb  and  70 
mb.  Figure  B3  shows  the  cross-sections  that  are  taken  for  the  graphics  that  follow  in 
Sections  3. 2-3. 9.  Of  note,  the  only  cross-section  plotted  from  Figure  B3  is  the  cross- 
section  labeled  M2. 
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3.2.  Parameter  Extraction  and  Reynolds  Decomposition 

As  can  be  seen  from  Equation  (B8),  several  different  variables  need  to  be  extracted 
from  the  model  output  in  order  to  perform  a  preliminary  TKE  budget  analysis.  The 
STRATONHMASS  model  was  formatted  to  output  the  U-relative,  V-relative,  and  W- 
relative  wind  field  in  accordance  with  S88’s  restriction  that  the  coordinate  system  be 
aligned  with  the  mean  wind  field.  The  u,  v,  and  w  model-relative  winds  are  then  averaged 
over  the  3  5 -minute  simulation  using  Equation  (B2).  These  terms,  denoted  as  the  mean  u, 
V,  and  w  wind  fields,  are  used  in  Terms  1, 2,  3,  and  5  from  Equation  (B8). 

In  order  to  calculate  perturbations  from  this  3  5 -minute  average,  Equation  (B4)  was 
used  to  produce  a  perturbation  data  set,  denoted  u’,  v’,  and  w’,  from  the  respective  fields 
and  averages  calculated  using  Equation  (B2).  Figures  B4  and  B5  show  the  M2  cross- 
section  for  the  u’  and  w’  fields  at  forecast  time  t  =  19.  The  perturbation  fields,  when 
looped  for  each  model  simulation  time,  show  a  vertically  propagating  wave  coming  from 
below  300  millibars  and  moving  into  the  stratosphere.  The  full  physical  interpretation  of 
these  fields  is  deferred  to  future  research. 


Figure  B3.  Domain  map  showing  the  location  of  various  cross-sections  for  the  222-m 
STRATONHMASS  simulation  of  the  12  December  2002  case.  For  the  purposes  of  this 
paper,  the  only  cross-section  that  is  used  is  M2,  in  the  center  of  the  map. 
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Figure  B4.  U’  (color  shaded)  for  the  12  December  2002  simulation,  valid  at  2238  UTC, 
12  December.  The  main  stratospheric  wave  is  located  roughly  corresponding  to  the 
positive/negative  couplet  centered  between  1 50-80  millibars. 
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Figure  B5.  As  with  Figure  B4,  except  for  w’.  The  wave  tilt  is  clearly  seen  here  in  the 
vertical  velocity  perturbations,  with  a  very  nice  positive/negative/positive  couplet  phase 
tilted  with  the  stratospheric  wave. 


The  same  Reynolds  decomposition  is  applied  to  the  virtual  temperature  field  from 
the  MASS  model  output.  First,  the  virtual  temperature  field  is  converted  to  virtual 
potential  temperature,  and  then  decomposed  using  Equations  (B2)  and  (B3).  Figure  B6 
shows  the  virtual  potential  temperature  perturbation  (9v’)  for  the  12  December  case  at 
time  t  =  19,  the  same  time  as  shown  in  Figures  B3  and  B4.  The  large  9v’  corresponds  to 
the  same  location  where  the  stratospheric  wave  is  increasing  the  heights,  and  the  small  Bv’ 
corresponds  to  the  location  where  the  wave  is  decreasing  the  heights.  Further  research 
must  be  done  on  this  hypothesis,  but  quick  visual  inspection  and  conferencing  with  Dr. 
Michael  Kaplan  and  Paul  Suffern  aided  this  argument. 

3.3.  TKE  Tendency  Due  to  Terms  1  and  2:  U  and  V-Advection 

Using  Equation  (B8),  Term  1  and  Term  2  were  calculated  using  the  time-averaged 
u  and  V  mean  wind  field  and  GEMPAK’s  horizontal  derivative  scheme  (DDX  and  DDY) 
for  the  TKE  field.  Following  Equation  (B3),  the  average  TKE  field  was  also  calculated 
for  the  3 5 -minute  simulation.  Thus,  the  contributions  from  Term  1  and  Term  2: 
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can  be  quantified  over  the  35-minute  model  simulation. 

Figure  B7  shows  Term  1  plotted  for  the  M2  cross  section  averaged  over  the  350- 
minute  domain. 


Figure  B6.  As  with  Figures  B4  and  B5,  but  for  perturbation  virtual  potential  temperature. 
The  field  is  a  bit  higher  up  than  the  u’  and  w’  maxima,  but  is  of  the  same  shape  and 
orientation  at  this  time. 
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Figure  B7.  Term  1  (U-Advection)  plotted  for  the  3 5 -minute  time  averaged  field.  Notice 
the  lone  maxima  and  minima  couplet  centered  around  70  mb,  much  higher  than  the  actual 
stratospheric  wave. 


Term  1  is  relatively  small  everywhere  within  the  domain  expect  near  the  apex  of 
the  breaking  stratospheric  wave,  which  is  zoomed-in  and  shown  close-up  for  Figure  B8. 
These  results  are  consistent  with  the  findings  of  S88,  and  Term  1  is  likely  quantitatively 
small  compared  to  the  other  terms. 

Similar  findings  are  found  in  Term  2,  shown  in  Figure  B9.  Like  Term  1,  Term  2  is 
small  except  for  small  pockets  in  the  upper  troposphere  where  the  stratospheric  wave  has 
broken.  Term  2  is  also  likely  relatively  small  compared  to  the  other  terms  in  the  TKE 
budget.  Further  analysis,  particularly  a  small  time  average  from  Equation  (B3),  is  needed 
and  is  planned  in  future  research. 

3,4.  TKE  Tendency  Due  to  Term  3:  W-Advection 

Using  Equation  (B8),  Term  3  was  calculated  using  a  separate  Fortran  program, 
written  by  the  author.  The  program  calculates  the  vertical  derivative  of  turbulent  kinetic 
energy,  which  has  been  averaged  over  the  3 5 -minute  mode!  run.  From  Equation  (B9): 
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Using  the  w  average  wind  field,  Term  3  can  be  calculated  using  Equation  (B8); 
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8t  dp 

Figure  BIO  shows  the  35-minute  averaged  cross-section  for  Term  3.  Because  w  is  being 
used  here  instead  of  omega,  the  plot  is  strictly  qualitative.  The  first  set  of  large  Term  3 
values  are  shown  zoomed  in  Figure  B 11,  in  a  layer  between  300  and  250  millibars.  The 
area  of  strongest  vertical  motion  associated  with  the  convection  lies  just  below  this  layer 
and,  thus,  the  convection  itself  drives  higher  values  of  Term  3  here.  The  second  set  of 
large  Term  3  values  are  shown  zoomed  in  Figure  B12,  in  the  lower  stratosphere  between 
120  and  65  millibars.  Here,  the  breaking  stratospheric  wave  has  created  a  large  vertical 
gradient  in  turbulent  kinetic  energy. 


Figure  B8.  Term  1,  U-Advection.  As  with  Figure  B7,  but  zoomed  to  a  vertical  axis  of  50- 
100  mb. 
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Figure  B9.  Term  2:  V-Advection, 
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Figure  BIO.  Term  3:  W-Advection. 
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Figure  B1 1.  Term  3,  W-Advection,  with  a  vertical  axis  from  180-280  millibars. 
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Figure  B12.  Term  3,  W-Advectiorij  but  with  a  vertical  axis  from  60-110  millibars. 


Figure  B13  shows  the  same  color-filled  values  for  Term  3,  with  contoured  values 

dc 

of  — .  Clearly,  the  large  Term  3  values  are  driven  primarily  because  of  a  high  gradient  in 

yp 

dc 

—  rather  than  a  flux  of  w-momentum.  The  area  of  highest  Term  3  lie  above  the 

dp 

stratospheric  wave  break,  leading  the  author  to  hypothesize  that  the  rapid  change  from 
large  to  small  TKE  values  above  the  wave  (de  <  0  and  dp  <  0)  acts  to  generate  higher 
TKE  generation  there.  This  result  will  be  tested  in  future  research,  as  aircraft  who  choose 
to  fly  over  convection  may  actually  be  flying  into  areas  of  high  TKE  generation  in  the 
lower  stratosphere  holds  in  further  applications. 

3.5.  TKE  Tendency  Due  to  Term  4:  Buoyancy 

Using  Equation  (B8)  and  the  perturbation  fields,  Term  4  is  calculated  using  the  35- 
minute  average  virtual  potential  temperature,  and  the  product  of  w'  and  9v’.  From 
Equation  (B8); 
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—  »oc...X(w'^/) ...  TERM  A 

dt  e. 

The  quasi-heat  flux  w’0v’  was  subsequently  averaged  over  the  35-minute  model  run  to 
produce  Figure  B14,  which  is  Term  4  plotted  in  the  same  1VI2  cross-section  as  the  other 
figures.  Figure  B14  shows  an  interesting  couplet  pattern  centered  around  90  millibars; 
about  15-20  millibars  lower  than  the  highest  values  in  Term  3.  Further  investigation  of 
these  phenomena  will  be  required  to  gain  a  physical  understanding  to  the  couplet’s 
meaning. 


Figure  B13.  Term  3,  W-Advection,  with  contoured  values  of  de/dp. 
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Figure  B14.  Term  4,  Buoyancy  production. 


3.6.  TKE  Tendency  Due  to  Term  5:  Mechanical  Shear 

Using  Equation  (B8),  and  the  perturbation  fields,  Term  5  is  calculated  using  the  u’ 
and  w’  fields  multiplied  together  and  averaged  over  the  35-minute  model  run  using 
Equation  (B2).  The  average  basic  state  shear  of  the  u-mean  wind  is  calculated  using  a 
separate  Fortran  program  and  has  different  resolution  above  and  below  100  millibars  as 
described  earlier  in  this  paper.  From  Equation  (B8),  Term  5  is  expressed  as: 


—  «oc ...  +  (u'w')pg —  ...  TERMS 
■  dt  dp 

Figure  B15  shows  the  35-minute  average  for  Term  5  using  the  same  M2  cross-section  as 
the  other  figures. 

The  most  noteworthy  information  about  Figure  B 15  is  the  anomalously  high 
values  of  Term  5  throughout  a  large  part  of  the  domain.  At  the  present  time,  the  author 
speculates  that  the  coordinate  transformation  may  be  the  culprit,  however,  further 
investigation  is  necessary  to  be  certain.  However,  though  the  magnitude  is  much  too 
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large,  the  pattern  of  mechanical  shear  production  is  consistent  with  theory.  Since  the 
stratosphere  is  extremely  stable,  one  suspects  that  the  dominant  production  term  for  the 
TKE  tendency  equation  would  be  mechanical  shear,  due  to  the  large  gradient  in  wind 
speed  produced  by  the  wave. 

3.7.  TKE  Tendency  Due  to  Term  6;  Turbulent  Transport 

Using  Equation  (B8)  and  the  perturbation  fields,  Term  6  is  calculated  using  the  w’ 
and  turbulent  kinetic  energy  (e)  fields  multiplied  together  and  averaged  over  the  35- 
minute  model  run  using  Equation  (B2).  The  quantity  w’e  represents  the  turbulent  flux  of 
TKE  from  the  vertical  velocity  perturbation.  From  Equation  (B8),  Term  6  is  expressed  as: 


de 


dw'e 
+ - 


dp 


TERM6 


Figure  B16  shows  the  35-minute  average  for  Term  6  using  the  M2  cross-section. 

There  are  two  primary  regions  of  interest  in  Figure  B16.  First,  several  large  Term 
6  maxima  appear  near  100  mb.  As  speculated  in  Section  2,  Term  6  is  likely  largest  where 
gradients  in  large  values  of  w’  and  e  are  juxtaposed.  Further  investigation  of  the 
magnitude  of  w’  and  e  in  this  layer  are  necessary  to  determine  why  there  are  positioning 
maxima  there.  In  addition,  a  couplet  in  Term  6  appears  on  the  eastern  half  of  the  domain 
in  far  western  Louisiana.  According  to  discussion  with  Mr.  Suffem,  convection  is  ongoing 
in  this  region  not  tied  directly  to  the  vertically  propagating  gravity  wave. 
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Figure  B15.  Term  5,  Mechanical  shear  production. 
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Figure  B16.  Term  6,  Turbulent  transport. 


3.8.  TKE  Tendency  Due  to  Term  7:  Pressure  Perturbations 

Using  Equation  (B8)  and  the  perturbation  fields,  Term  7  is  calculated  using  the  w’ 
and  p’  fields  multiplied  together  and  averaged  over  the  35-minute  model  run  using  (2.2). 
The  pressure  perturbation  term  is  outputted  directly  from  the  STRATONHMASS 
simulation.  The  quantity  w’p’  represents  the  creation  and  destruction  of  TKE  from  the 
pressure  perturbations.  From  Equation  (B8),  Term  7  is  expressed  as: 


de 

dt 


fBOC...  +g 


dw'p' 


TERM? 


Figure  B17  shows  the  3  5 -minute  average  for  Term  7  using  the  M2  cross-section.  The 
blotchy  nature  of  the  cross-section  is  partly  to  blame  for  the  appearance  of  the  data. 
However,  there  appear  to  be  several  couplets  in  Term  7  located  in  the  same  location  as 
other  Terms  in  the  TKE  equation.  As  with  the  other  terms,  further  analysis  is  needed  to 
explain  the  pattern. 


66 


Figure  B17.  Term  7,  Pressure  perturbations. 


3.9.  TKE  Tendency  Due  to  Term  8:  Dissipation 

Using  the  same  Reynolds  decomposition  algorithm  described  in  sections  2. 1-2.3, 
the  average  turbulent  kinetic  energy  dissipation  was  averaged  over  the  3  5 -minute  model 
period.  TKE  is  a  dissipative  quantity  and,  hence,  in  the  absence  of  Terms  1-7,  TKE 
always  acts  to  be  reduced.  From  S88,  TKE  dissipation  is  highest  in  areas  where  TKE  is 
large,  since  overall,  TKE  is  conserved.  Figure  B18  shows  the  3  5 -minute  average 
dissipation  from  the  12  December  case. 

Term  8  is  surprisingly  small  in  relation  to  the  other  terms.  This  suggests  that  the 
explicit  calculation  of  the  terms  in  the  TKE  budget  may  be  anomalously  high  on  many 
different  levels.  The  possibility  also  exists  that  the  model  doesn’t  resolve  dissipation  as 
well  away  from  the  PEL  as  is  expected  by  calculating  this  explicitly. 
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Figure  B18.  Term  8,  TKE  dissipation. 


4.  FUTURE  WORK 

As  highlighted  in  many  sections  of  this  paper,  there  is  much  work  to  be  done  on 
the  TKE  budget  algorithm  and  output  for  the  cases  prescribed  in  the  Air  Force  project. 
The  following  research  objectives  for  the  Summer  2005  term  are  proposed; 

•  Modify  current  existing  model  output  for  height  coordinates,  and  utilize  Equation 
(B5),  the  TKE  tendency  equation  in  height  coordinates.  This  requires  that  the 
MASS  model’s  vertical  resolution  for  the  stratospheric  data  are  fixed,  and  set 
intervals  are  prescribed  in  order  for  the  vertical  derivatives  to  be  consistent; 

OR 

•  Re-compile  the  u  and  v  relative  wind  fields  in  such  a  fashion  that  the  GEMPAK 
program  GDOMEG  can  be  used  to  calculate  an  approximate  omega  vertical 
velocity  using  the  continuity  equation.  The  modified  TKE  tendency  in  pressure 
coordinates,  given  in  Equation  (B8},  can  then  be  used  for  complete  quantitative 
investigation.  A  further  analysis  of  the  approximations  that  are  invoked  will  be 
required.  This  path  may  be  able  to  give  a  quicker  mathematical  interpretation  than 
the  previous  suggestion; 
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•  Increase  the  order  of  the  finite  difference  scheme  employed  to  calculate  the 
vertical  derivatives  and  compare  to  the  first-order  scheme; 

•  A  full  quality  assurance  check  of  all  the  data  ingested  by  the  TKE  budget  code  to 
assure  that  no  bad  data  are  responsible  for  strange  budget  calculations  (see  Section 
3.9); 

•  Compare  the  TKE  tendency  equation  results  to  that  of  the  model’s  TKE  average 
over  the  averaged  time  period  (this  will  require  a  direct  computational  TKE 
budget).  The  author  and  colleagues  will  then  investigate  the  role  of  each  term  in 
generating  TKE  tendency  physically  using  wave  dynamics  and  previous  literature; 

•  Complete  the  TKE  budget  for  all  Air  Force  project  cases  and  regimes  as  prescribed 
by  Dr.  Michael  Kaplan; 

•  Use  the  TKE  budget  analysis  code  and  script  to  develop  an  EDR  budget  algorithm 
that  closely  resembles  that  of  the  TKE  budget.  The  EDR  tendency  equation  (not 
shown)  uses  many  of  the  same  variables  as  the  TKE  budget  and  should  not  require 
a  great  deal  more  in  the  way  of  calculations  and  manipulations.  The  EDR  budget 
would  then  be  calculated  for  all  the  Air  Force  project  case  studies. 

5.  CONCLUSIONS 

A  preliminary  TKE  budget  analysis  has  been  presented  for  the  12  December  2002 
large-amplitude  gravity  wave  case  for  a  model  simulation  time  period  of  35  minutes  at  a 
222-m  horizontal  resolution.  Due  to  various  data  problems,  the  TKE  tendency  equation 
was  written  in  pressure  coordinates  and  rather  than  height  coordinates.  The  coordinate 
transformation  necessary  to  perform  many  of  the  calculations  leaves  the  budget  to  become 
strictly  a  qualitative  one. 

The  purpose  of  this  project  was  to  develop  the  code  and  perform  an  initial  analysis 
of  the  TKE  budget  for  one  case  study  as  a  stepping-stone  to  further  research.  The  results 
presented  in  Section  3  are  purely  of  a  qualitative  nature  due  to  the  various  approximations, 
data  ingestion  issues,  and  coordinate  transformations.  The  development  of  a  skeleton  code 
in  which  TKE  budgets  can  be  calculated  is  of  great  significance  to  the  project  and  future 
research  described  in  Section  4. 

The  TKE  budget  analysis  raises  many  research  questions  that  will  require  quite  a 
bit  of  thinking.  Are  the  patterns  presented  in  the  terms  of  physical  significance  or  strictly  a 
byproduct  of  the  calculations  within  the  budget?  Many  more  TKE  budgets  must  be 
performed  on  the  various  case  studies  in  order  to  begin  answering  some  of  these 
questions.  A  rigorous  physical  interpretation  of  the  results  is  deferred  from  this  conclusion 
section  to  further  research. 

6.  SUMMER  2005  RESEARCH  UPDATE 

As  discussed  in  Section  4,  Future  Work,  several  different  research  goals  were 
described  for  the  Summer  2005  term.  The  first  objective  involved  utilizing  (2.4)  to 
transform  the  TKE  budget  into  a  height  coordinate  form.  After  the  TKE  budget  in  pressure 
coordinates  was  examined  closer,  it  was  confirmed  that  the  units  were  not  consistent  in  the 


69 


derivation.  The  value  of  the  pressure  coordinate  exercise  was  to  develop  the  first  skeleton 
code  required  for  the  budget  and  identify  the  regions  of  interest. 

The  TKE  budget  code  was  then  developed  for  use  in  a  height  coordinate  system 
consistent  with  the  PBL-inspired  Equation  (B4).  Included  in  this  Appendix  is  the  71-m 
and  222-m,  12  December  2002,  TKE  budgets  in  height  coordinates.  Cursory  examination 
of  the  71"  and  222-m  TKE  budgets  indicate  a  finer-scale  resolved  turbulence  signal  within 
the  71-m  run,  as  well  as  stronger  emphasis  on  the  u  and  v  advection  terms.  The  increase  in 
TKE  contribution  from  the  u  and  v  advection  terms  indicate  a  stronger  horizontal  TKE 
gradient  within  the  71-m  run,  consistent  with  the  finer  horizontal  resolution  and  the 
expected  resolve  turbulence  signal.  The  term  calculations  for  the  71-  and  222-m  runs  are 
shown  in  Figures  B19-B34.  Future  analysis  and  discussion  is  deferred  to  future 
documents  due  to  time  constraints. 
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Figure  TKE  Budget  for  12  December  2002, 222-m  case  in  height  coordinates: 

Term  1,  U-Advection. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B20.*  TKE  Budget  for  12  December  2002, 222-m  case  in  height  coordinates: 
Term  2,  V-Advection. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B21.*  TKE  Budget  for  12  December  2002, 222-m  case  in  height  coordinates: 
Term  3,  W-Advection. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B22.*  TKE  Budget  for  12  December  2002, 222-m  case  in  height  coordinates: 
Term  4,  Buoyancy  Production. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B23.*  TKE  Budget  for  12  December  2002, 222-m  case  in  height  coordinates: 
Term  5,  Mechanical  Shear. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B24.*  TKE  Budget  for  12  December  2002, 222-m  case  in  height  coordinates: 
Term  6,  Turbulent  Transport. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 


75 


21000  - 


20000  - 


19000  - 


18000  - 


17000  - 


16000  ' 


15000  - 


14000  - 


13000  ^ 


12000  - 


11000  - 


lOOOO- 


9000— p-r- 
30.40;-9i55 


-O.flCl  -7.5i-*^E.ti-C<-S.Si-fl4  t.eCO  E.5E-C4  5.Cs-(14  V.5e-I)4  0.001 


Figure  B25.*  TKE  Budget  for  12  December  2002,  222-m  case  in  height  coordinates; 
Term  7,  Pressure  Perturbations. 


^  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update, 
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Figure  B26.*  TKE  Budget  for  12  December  2002, 222-m  case  in  height  coordinates: 
Term  8,  Model  Derived  EDR. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B27.*  TKE  Budget  for  12  December  2002, 71-m  case  in  height  coordinates: 
Term  1,  U-Advection. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B28.’^  TKE  Budget  for  12  December  2002, 7i-m  case  in  height  coordinates: 
Terra  2,  V-Advection. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B29.*  TKE  Budget  for  12  December  2002,  71-m  case  in  height  coordinates: 
Term  3,  W-Advection. 


Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B30.*  TKE  Budget  for  12  December  2002,  71-m  case  in  height  coordinates: 
Term  4,  Buoyant  Production. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B31.*  TKE  Budget  for  12  December  2002,  71-m  case  in  height  coordinates: 
Term  5,  Mechanical  Shear. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B32.*  TKE  Budget  for  12  December  2002,  71-m  case  in  height  coordinates: 
Term  6,  Turbulent  Transport. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B33.*  TKE  Budget  for  12  December  2002, 71-m  case  in  height  coordinates: 
Term  7,  Pressure  Perturbations. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 
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Figure  B34.*  TKE  Budget  for  12  December  2002, 71-m  case  in  height  coordinates: 
Term  8,  Model-Derived  EDR. 


*  Denotes  figures  inserted  as  part  of  Summer  2005  Research  Update. 


85 


REFERENCES 


Bluestein,  Howard  B.,  1993:  Synoptic-Dynamic  Meteorology  in  Midlatitudes,  Volume 
II.  Oxford  University  Press:  Chap.  3.4,  pp.  431-520. 

Lin,  Y-L.  Mesoscale  Meteorology.  To  be  published.  Ch.4:  Finite  Difference  Schemes, 
Spring  2005. 

Stull,  Roland  B.,  1988.  An  Introduction  to  Boundary  Layer  Meteorology.  Khmer 
Academic  Press:  639  pages. 


87 


